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ABSTRACT 

Electron energy distribution functions have been calculated in a 
U -plasma at 1 atmosphere for various plasma temperatures (5000-8000“K) 
and neutron fluxes (2 x lQ^^-2 x 10^^ neutrons/ (cm^-sec)). The distribu- 
tions are assumed to be a summation of a high-energy tail and a Maxwel- 
lian distribution. The sources of energetic electrons considered are the 
fission-fragment induced ionization of uranium and the electron induced 
ionization of uranivim. The calculation of the high-energy tail is reduced 
to an electron slowing down calculation, from the most energetic source 
(~ 2.1 keV) to the energy where the electron is assianed to be incorporated 
into the Maxwellian distribution (~ 15 eV) . The pertinent collisional 
processes are electron-electron scattering and electron induced ionization 
and excitation of uranium. 

Two distinct methods have been employed in the calculation of the 
distributions. One method is based upon the assumption of continuous 
slowing and yields a distribution inversely proportional to the stopping 
power. An iteration scheme is utilized to include the secondary electron 
avalanche. 

In the other method, a governing equation is derived without assuming 
continuous electron slowing. This equation is solved by a Monte Carlo 
technique which simulates Coulombic collisional slowing analytically while 
ionization and excitation events are simulated in a random \ifalk fashion. 
Consequently, the secondary electron avalanche is included explicitly. 

Both methods yield comparable results at high energies 100 eV) , 
with disparities arising at lower energies due to the inapplicability of 
the continuous slowing assumption. The distribution functions calculated 
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in both models are observed to be linearly dependent upon the neutron flux 
while inversely proportional to the plasma temperature. The electrons 
within the calculated high-energy tail induce ~ 10^^ more excitations of 

3 

uranium per cm per second than are induced by Maxwellian electrons. 

Since the threshold of non-Maxwe Ilian behavior is ~ 15 eV, the present re- 
sults suggest seeding the plasma with a species having a high excitation 
threshold, e.g. helium, in order to better capitalize upon the excitation • 
characteristics of the high-energy tail in possible applications as a 
lasing medium or a radiation source. 
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CHAPTER I 
INTRODUCTION 


A. Definition of Problem 

The objective of this study is the. deduction of the effect of the 
presence of fission-fragments upon the electron energy distribution 
function in a uranium plasma. Parametric studies of the neutron-flux 
(a measure of the fission-fragment density) and temperature dependence 
of the distribution function are xindertaken to provide insight into the 
plasma conditions- under which the fission-fragments have the most pro- 
nounced effect upon the distribution function. These calculations will 
be obtained by two separate models; a simple model for survey calcula- 
tions and a second, more refined model, by which the accuracy of the 
former may be judged. 

Primary emphasis is placed upon the high-energy tail of the dis- 
tribution function, i.e., at energies above the excitation threshold, as 
it is anticipated that the results generated herein will be used for the 
calculation of excitation rates. At such energies, the calculation of 
the distribution function reduces to the problem of slowing down from a 
source. The source of electrons to be considered here is comprised of 
two distinct components, each distributed in energy. The first consists 
of those electrons generated by the fission-fragment induced ionization 
of the background uranium during their slowing down process, while the 
second consists of those secondary electrons produced through the ioniza- 
tion of uranium by energetic electrons as they thermaliz'e. 
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B. Motivation 

The motivating forces behind this study can best be identified by 
examining some of the anticipated applications of a uranium plasma. His- 
torically, the first application envisioned was the- utilization of the 

plasma as an energy source for an ion rocket engine termed’ the nuclear 

fl-4‘) 

light -bulb concept.*- . The success of the idea depends upon the ef- 
ficiency at which the energy released during the fissioning of uranium 
is transmitted to and absorbed by the hydrogen fuel.. The energy trans-.. 
mission process consists of a conduction chain and a radiation chain. 

In the conduction chain, the energy deposited within the plasma by the 
neutrons and fission-fragments is conducted away from its source to the 
hydrogen fuel. In the radiation chain, a portion of the fission-fragment 
energy is transmitted to electrons through the ionization of the background 
uranium. The electrons in turn excite the background which transforms 
the energy into radiation as the atoms de-excite. Then both the line 
radiation and the blackbody radiation pass through a ’’window" into the 
fuel where it is to be absorbed. 

Two additional applications make extensive use of the radiation 
chain, namely direct nuclear pumping^^"^^ and photo-chemical production 
by extracting energy from the uranium plasma in the form of light. In 
the scheme of direct nuclear pumping, a population inversion is sought 
either by seeking a situation where the uranium will lase or by transfer- 
ring the energy from the uranium plasma to a second species which would lase. 
In the latter scheme, less stringent requirements are placed upon the exci- 
tation rates, as a population inversion of the uranium itself is unnecessary. 
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The success of this scheme dep'ends solely upon the ability to shape the 
radiation spectrum via either the determination of the plasma opacity or 
the- transluscent properties of the "window". 

A fourth application involves the high efficiency extraction of 
energy from the uranium plasma through an MHD cycle. This scheme may 
be combined with either of the two previous schemes which would serve as 
topping cycles to further enhance efficiency. 

From these possible applications, the importance of the calculation 
of the electron energy distribution function can be gauged. The success 
of most of these schemes depends upon an accurate detemiination of the 
excitation rates. A prime means of exciting atoms is through electron 
induced excitation. Consequently, at the heart of the problem is the 
need for a detailed knowledge of the number of electrons capable of in- 
ducing excitation. Such a query can be satisfied only with a detailed 
calculation of the electron energy distribution function rather than 
assuming erroneously that the distribution is Maxwellian in radiation 
calculations . 

C. Description of Plasma 

1. Classification of Plasma 

The plasma conditions to be investigated include temperatures 

ranging from 5000°K to 8000°K (the boiling point of uranium is' 4407°K 

at one atmosphere); a pressure of one atmosphere, and neutron fluxes 

ranging from 10 to neutrons/ (cra^-sec) . In determining the rate .of 

occurrence of fission reactions,, the uranium is assumed to consist en- 
235 

tirely of the U isotope. Furthermore, the neutrons are assumed to be 
in thermal equilibrium with the plasma so that the fission cross-section 



4 


f33') 

calculated by Bussard'- *' is applicable, i.e., the effects of neutron 
spectrum hardening are negligible. In his calculation, the fission 
cross section weighted by the neutron flux distribution is averaged over 
energy, thereby eliminating the energy dependence of the cross section, in 
favor of a characteristic temperature. 

Upon fissioning, a uranium atom is assumed to split into two fission- 
fragments. The lighter fragment (96 amu) is bom at 98. MeV and a charge 
of +16e, while the other is bom at 67 MeV and a charge of :^■15e. The 
distribution of each fission- fragment is taken .to be inversely proportional 
to energy. The tendency of a fission-fragment to neutralize its 
positive charge as it thermalizes is included by assuming the fission- 
fragment's charge to be proportional to velocity. 

In terms of the previously described energy extraction schemes, the 
plasma conditions cited above are characteristic of a subcritical uranium 
plasma. Also, this temperature range encompasses the 6000®K temperature 
of a proposed 5-MlV experimental reactor, The conditions of ,a critical 
plasma are somewhat hotter (center line temperature of 4,0,000"K^^^^i) and 
of higher pressure (approaching 500 atmospheres) . However, the results 
for the plasma conditions to be studied should be applicable to the outer 
boundary layer of a critical plasma. 

For the plasma parameters cited, the densities of the various plasma 

constituents can be predicted by the Saha equations provided the 

necessary partition functions are known. Due to the lack of experimental 

data, the ratio of the partition functions is assumed to be unity (after 
( 12 ) 

Krascella ). The results for the Saha predicted densities appear in 
Fig. 1. A first order approximation of the perturbation to these densities 



TEMPERATURE, °K 


Fig. 1. The density of the uranium plasma constituents plotted versus 
temperature assuming a pressure of one atmosphere. 
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caused .by the production o£ fission-fragment generated electrons is only 
of the order of 1/10% for cases of interest here. Therefore, a further 
correction for radiation effects is generally negligible. 

•With the electron densities and temperatures incurred, the plasma 
defies classification in classical terms (see Fig. 2) . There exists an 
insufficient number of charged particles within a Debye sphere to provide 
the necessary screening of a charged test-particle required for the sat-^ 
isfaction of the binary collision ‘assumption of the classical hinetic plas- 
ma. Similarly, the Landau distance of correlation is not sufficiently 
large for the plasma to be characterized by the very strong correlations 
of the classical collective plasma. Due to these difficulties, both a 
binary collisional treatment and a unified treatment which in- 
corporates collective interactions as well as the binary collisions (see 
Appendix A) are applied to the Coulombic collisions. 

2. Delineation of Collisional Processes 

The dominant types of electron collisions present in the uranium 
plasma are: the aforementioned elastic Coulombic collisions — electron- 

electron and electron-ion scattering; and the inelastic collisions — , 
ionization and excitation of neutral and singly ionized uranium (see 
Appendix B). 


Since the inelastic cross sections. have not been measured experi- 
mentally, they must be calculated from formulae based upon the Gryzinski 
model, implementing the ionization and excitation data of Parks, et 

al.^175 



Fig. 2. A delineation of the plasma parameters defining the classical 
kinetic plasma and the classical collective plasma. The Debye 
length is given by h, the electron density by n^, the charged 
particle spacing by dg and the Landau correlation distance by r^ 
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D. Methods o£ Solution 

1 . Governing Equation 

The establishment of a governing equation for the distribution 
function is complicated by' the presence of both inelastic and Coulombic 
collision types. For plasmas dominated by just one of these types, there 
exists a plasma equation describing the collision mechanisms involved. 

The Boltzmann equation, based upon the assumption that the duration of 
a collision is much less than the time between collisions, depicts binary 
collisions, short range forces, and neutral scattering; whereas the 
Fokker-Planck equation with a cut-off distance depicts Coulombic col- 
lisions which are not strictly binary in nature due to the long-range 
force responsible for the Coulombic interaction. 

' Then, the governing equation must be a' combination of these two 

equation types in order to accommodate the presence of both collision 

ri91 

types. Such an equation has been formulated by Dreicer for a partially 
ionized gas. However, he included electron-neutral collisions which are 
negligible in the present case. 

The distribution functions derived from both of these equations 
either separately or combined are one particle distributions. Such dis- 
tributions fail to describe the correlation effects anticipated in a 
uranium plasma. They can only .be described accurately by a many-bodied 
distribution which satisfies the Liouville equation. The complexity of 
this latter equation renders it impractical for direct use. However, 

Upon employment of the Coulombic cross section, the Boltzmann equation 
reduces to a Fokker-Planck type equation which is equivalent to the 
Fokker-Planck equation under special circumstances (see Montgomery and 
Tidman C20) ^ _ 
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the correlation effects can be approximately described within a specified 
error tolerance via the BBGKY hierarchy. The Boltzmann and Fokker- 
Planck equations represent a zero order kinetic equation in this hierarchy 
while the, first orHer equatioh’ l's' ‘the Lenard-Balescu equation.. .^The 
accuracy of these equations is expressed in terms of the assumed small 
plasma parameter g = ^ ^ « 1, where n^ is the charged particle density 

and Ajj is the Debye length. The zero order kinetic equations are ac- 
curate to order 1, the first order kinetic equation is accurate to order 
g, the second order equation is accurate to order g"", etc. In the uranium 
plasma, the assumption of g being small is not well satisfied. This pre- 
sents a problem since the accuracy of even higher order equations becomes 
uncertain. Consequently, two sets of governing equations will be separ- 
ately imposed; they are a combination Boltzmann and Pokker- Planck equa- 
tion and a combination Boltzmann and Lenard^-Balescu equation. It is 
argued that a comparison of the resuJts wi.il provide an estimate of the 
error introduced by not employing a many-bodied distribution function. 

The previous equations may he simplified by noting tlie uranium 
X^lasma to be in a steady state, Tiic presence of Lite higii- energy electrons 
produced botli by f Is si on- fragments and other high-energy electrons creates 
a non -equilibrium state. However, the results of and al 

similar electron source rates indicate tho source electrons relax into a 
Maxwellian distribution, with the non- equilibrium effects restricted to 
high energies. Then, the problem becomes one of investigating the re- 
laxation of the high-energy tail into a Maxwellian distribution as de- 
scribed by the collision terms of the aforementioned equations. 
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2. Review o£ Methods and Proposed Solution 

Most methods appearing in the literature are applicable to one 

equation type. Most notable of the methods are those of Rosenbluth, 

f241 r25') 

et al.^ ^ for the Fokker-Planck equation and those of Nighen^ ^ and 

Holstein for the Boltzmann equation. They all ejq)and the distribution 
function in terms of Legendre polynomials. Such methods are ideally 
suited to anisotropic plasmas with E-fields or injected beams. 

in the present case, the assumption of kn. isotropic source and a 
primary interest in the high-energy tail make the method of Fano^^^’^^^ 
much more appealing. However, provisions must be made to include a nascent 
or fission-fragment generated source distributed in energy plus a .secondary 
electron source. The tractability of the resulting solution for the dis- 
tribution fimction renders it an -important tool for survey calculations 
and in the analysis of the distribution function. 

In order to partially relax the assumption of continuous slowing 
down inherent in Fano's method, a Monte Carlo simulation is also per- 

foimied which separates the Coulombic collisions from the- inelastic col- 

(29') 

lisions in a manner analogous to that of Wells ^ ^ in earlier analytic 

studies. Unlike the Monte Carlo calculation of Thomas and Thomas, 
the variation of the mean free path length with energy between the point 
of origin and the collision point is included in the calculation of the 
distance of random walk. Due to the presence of a secondary electron 
source, the ergodic hypothesis is' not applicable so that a number of 
particles correlated in time must be considered simultaneously rather 
than repeatedly simulating an individual electron. The increased degree 
of sophistication of this calculation is obtained at the ejq)ense of the 
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economy of the solution. Then, the practical application of this tech- 
nique would require that it be used solely as a check upon the validity 
of the analytic solution. 
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CHAPTER II 
ANALYTIC SOLUTION 


A. Introduction 

In this chapter, a simple analytic formulation is sought for the 
relaxation of supertherraal electrons into a thermalized ensemble of elec- 
trons which can be described by a Maxwellian distribution. The approach 
employed follows that of Spencer and Fano^^^^ for energetic electrons in 
an infinite medium, which predicts a distribution proportional to the 
inverse of the stopping power. Then, the resulting distribution will be 
a superposition of a high-energy tail on a Maxwellian thermal distribu- 
tion. From this formulation, the effect of varying several plasma para- 
meters can readily be predicted, or conversely, variations in the dis- 
tribution can easily be traced to their source. The ease of analysis 
afforded by this method renders it an important tool for surveys, but the 
more detailed calculation of Chapter III must be retained if high accuracy 
is desired. 

B . Derivation* 

An analytic solution for the electron energy distribution function 
in an infinite medium can be derived from the following, completely 
general, expression for the conservation of electrons in energy space: 

SCE)d£ + 

The distribution function f(E) represents the number of electrons per unit 
spatial volume per unit energy. The density of the electrons in the high- 

The derivation given here is a modification of one by Safanov. 
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energy tail can be obtained as follows: 



CE) JLe: = n 


C2) 


E ^ E-p 


where is the threshold for non-Maxwellian behavior of the distribution 
function. The source function S(E) is equal to the electron production 
rate from all sources per unit spatial volume per unit energy. The rate 
at which electrons recombine is represented by The rate of scat- 

tering from and into the energy interval dE about E is represented by 

and q^^, respectively. For the energies of interest, the electrons 
which are scattered into the energy interval dE about E originate at 
energies greater than E; that is, the upscattering of electrons can be 
assumed to be negligible. Then, the electron balance expressed in Eq. (1) 
can be visualized as in Fig. 3.- 

The scattering terms appearing in Eq. Cl) may now be evaluated. An 
e;q>ression for q^^^^ may be obtained by dividing f(E)dE, the number of 
electrons in the energy interval dE about E, by a characteristic deceler- 
ation time tg, i.e.. 


fa 


(3) 


rout "Tp" 

With the assumption of an infinitesimal energy loss per electron collision. 


dE/tg may be replaced by the rate at which electrons collisionally lose 
d.£ 

energy This is equivalent to the assumption of continuous slowing 
down. An additional implication of this assumption is that electrons 
scattered from one infinitesimal energy element dE about E + dE must 
be scattered into the adjacent element dE about E', i.e., 
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S(E)dE 



z/c(E)f(E)dE 


Fig. 3. The particle balance of Eq. (1) imposed upon the energy 
interval dE. about E. 
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fen = 

3 

where denotes the number of electrons/cm slowing into the energy in- 

terval dE centered about E, Substitution of Eqs. (3) and (4) into the 
balance of Eq. fl) gives: 


S(E)<lE- -h'P(E+c8£) =^-Fce)^! + 

i£+'4e (Atle 

or 

5 CE.) - CE^)f (.E) = ~ if CEf dLE.) 4E 1 - f (E) dEU 

^ dtlE^-a-E. dltk 

Invoking the fundamental theorem of calculus ^ Eq. (6) becomes 


S<.£) - = - 4 ^ f-fcEyisl \ 

4E AtleJ 


Integration of Eq. (7) yields 


^ OO 

J [ S<.E:) - ^(£') -P(E') j <!£.' = - (£' ) ASj \ 

In the limit as E approaches infinity, the distribution function f(E) 
vanishes. Consequently, the distribution function must obey the following 


Fredholm integral equation 


^Oo 

f(E)= ) [$(£')- 3^ (£')fcE')}4E.'/iE f- 

E ■’ dt £ 

At high energies (at or above the inelastic collision threshold), re- 


combination can be neglected. Defining the "cut-off energy" such that 
recombination ^d upscattering are negligible for E > E.^,, we obtain 
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■f(e) = 

Strictly speaking, Eq. (lOJ is valid only for the steady state of an 
infinite, isotropic medium where continuous slowing down is applicable 
and external forces or fields are absent. The assumption of continuous 
slowing down has proved to be a valid assumption for elastic scattering 
off of heavy targets xvhere the ratio of the energy lost to the original 
energy is small, as is evidenced by the successful application of the Fermi 
age theory^ '' to the slowing of neutrons by heavy moderators. This 
assumption should also be valid for the Coulombic collisions where small 
angle scattering (hence, small energy transfer) is dominant. In the 
present case, however, inelastic scattering, i.e., ionization and exci- 
tation collisions, represents an equally important energy loss mechanism 
v/hich does not necessarily comply with the continuous slowing assumption. 
This introduces some error into the model, thus, it can only be viewed as 
a first approximation. ‘ A more rigorous but more costly Monte Carlo 
ti-eatment is then developed in Chapter CIJ for more precise studies. 



C. Numerical .Solution 



Nujuerical results for the discnbution function can readily be ob- 
tained if ^ and S(E) are known. Provided the collisions are of a binary 
cIH 

nature, ^ can be decomposed into a sum of energy loss rates for each 


type of collision, i.e.. 


at at 


ORiQUfAL PAGj? tq 

*OOB Quattrr 


j Coulombic | ionizatioT 

collisions collision 


collision 
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Both ionization and excitation collisions qualify as binary collisions. 

However, Coulonibic collisions are not truly binary, yet they have been 

treated successfully as binary collisions in both applications of the 

f 191 

Fokker-Planck equation^ and in derivations of Fokker-Planck type 
equations. Then, the decomposition of the energy loss rate in 
Eq. (11) is applicable to the present case. 


dE 

Numerous treatments for jt- 

dt 

(14) 


as the Fokker-Planck model* 
and excitation may be obtained by 


= <£ >. 

d-t 


„ , , . exist in the literature, such 

Coulombic 

collisions 

The energy loss rates for ionization 


' I 


C12) 


where v is the speed of a test electron relative to thermal electrons 
[v corresponds to the energy E in Eq. (10)] and Z is the macroscopic 
inelastic cross section. The average energy loss per collision 


is defined as 




dE 


(13) 


C ) ,) rr 

where E is the energy of a test particle and E’ is the energy lost by 
the test particle as a result of a collision. The microscopic cross 
section o(E) and the energy transfer differential cross section — - ^ 

for excitation and ionization events necessary for evaluating Eq. (13) 
have not been heretofore measured experimentally nor calculated. Then, 
these quantities had to calculated specifically for this study from a 
Gryzinski model using the data of Parks, et al.^^^^ for uranium atom 


states (see Appendix B) , 
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2. Iterative Solution of Source Term 

The source term S(E) must include both secondary electrons and the 
nascent electrons resulting from ionization of uranium by fission-fragments. 
(Thermal electrons up-scattered in energy are excluded as they are neg- 
ligible for most of the energy interval of interest.) Since the dis- 
tribution of secondary electrons is dependent upon the distribution of 
nascent electrons and the manner in which they thermalize, a total source 
term S(E) cannot be known a priori. Consequently, S(E) and the distribu- 
tion function f (E) must be calculated in an- iterative manner as is de- 
scribed below. 

First, the production rate of nascent electrons designated by S^(E) 
is calculated (See Section D). These electrons relax into a primary 
electron distribution £q(E) according to Eq. (10). c During the thermali- 
zation process, the primary electrons further ionize the background 
uranium generating a source of secondary electrons Sj^(E). These secon- 
daries distribute themselves in energy as prescribed by Eq. (10), i.e., 
insertion of Sj^(E) in the equation yields f^(E), producing yet another 
generation of secondary electrons 82 (E) . This process is continued un- 
til the sum of the S^(E)'s converge to S(E) and likewise, the sum of 
the f ^ (E) ' s converge to f (E) . The convergence of the sum of the f ^ (E) ' s 

is readily obtainable within a few iterations, in agreement with earlier 

(32) 

observations of such a process by Fano and Spencer. 

D. Sample Results 

1 . Nascent Source 

The starting point in the iterative scheme to deteomiine the dis- 
tribution function is the Calculation of the nascent electron source S . 
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Since the electrons comprising are the result of fission-fragment 
induced ionization of uranium, it is essential that the fission-fragment 
distribution is knovm. A simple estimate of this distribution can also 
be had from Eq. (10) • The source of fission-fragments is so narrow in 
energy, it can be considered a delta function, i.e., two distinct 
fission-fragments are bom, as the result of a single fission event, at 
energies of 67 MeV and 98 MeV and masses of 140 amu and 96 amu, respectiv- 
ely. Then, Eq. (10) becomes 


f(E) = 

FF 


<lt 


(14) 


where S’ represents the number of fission events/ (cm^-sec) . Assuming a 

14 2 

neutron flux of 2 x 10 neutrons/ (cm -sec), an averaged fission cross 

(331 17 -3 

section of 57.6 bams, ^ and a gaseous uranium density of 5.6 x 10 cm 

q 

at 8000 K, S’ is evaluated to be 6.5 x 10 fission-fragments of each kind 

3 

are bom/ (cm -sec) in this example. 

The fission-fragments experience' electron capture over their entire 
track, i.e., q = < 1 qV/V^ where q^ and represent the initial charge 
(M6e) and velocity, respectively. Consequently, the energy loss dE/dx 
is a maximum at the beginning of their track. A semi-empirical formiila^^^^ 
for the energy loss of a fission-fragment at energy E is given by 

(15) 




da 


XC£.) 



where E^ is the energy at which the fission-fragment is bom and X is its 
range (see Eq. 3.50 of reference 10 for a semi-empirical expression for 
X) . Then, the fission- fragment distribution is 
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£(e) = 


XC£o) S- /fvj 


ZE. 


IE. 


where M is the mass o£ the fission-fragment and the relationship 
dE/dt = V dE/dx = • dE/dx is utilized. The nascent electron source 

appearing in Fig. 4 is then obtained by averaging the fission- fragment 
distribution over a Gryzinski energy transfer cross-section for ioniza- 
tion events, generalized for heavy, multi-charged ions.^^^^ The average 
energy of the nascent electrons is found to be '^^10 eV.. 


2 . Distribution Function 


The results of successive iterations upon the distribution function 
are also shown in Fig. 4. Convergence is easily realized in three iter- 
ations. The final solution of the high-energy tail (dot-dash Tine) -is ' 
displayed along with a Maxwellian distribution (solid line) corresponding 
to the plasma density and temperature previously cited. Where the high- 
energy tail intersects the Maxwellian, the source of electrons is no 
longer dominated by the nascent electrons and the'ir resulting avalanche 
but rather by up-scattered electrons. Therefore, it is assumed in Fig. ’4 
that the actual distribution will more likely resemble a summation of the 
Maxwellian and the high energy tail. 

3. Energy Loss Rate 

The energy loss rates necessary for the calculation of the distribu- 
tion function via Eq. (10) are displayed in Table 1. It is evident from 
the individual energy loss -rates that the Coulombic collisions are as im- 
portant in slowing down as are the inelastic- collisions, ionization and 
excitation, inspite of the' vast difference in the average energy lost per 
collision. Although the energy loss per collision by Coulombic interactions 
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Fig. 4. The nascent source versus energy as well as the distribu- 

tion function plotted versus energy. The distribution consists 
of a Maxwellian (“) and the' converged solution of the high-energy 
tail (-•-). The initial term in the series representation of fCE) 
is also plotted C — ) • 
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is small, the cross section is relatively large so that the loss rate 
for these events is comparable to inelastic loss rates. 

For high energies, the largest permissable energy loss per collision 
is a half of an electron’s energy before a collision. The probability 
of ■ such a hard collision, either Coulombic ox ionization, is small. 
Therefore, the assumption of continuous slowing is reasonable. However, 
at lower energies the energy lost in an ionization event can be a 
sizable fraction of the electrons original energy, and the continuous, 
slowing approximation becomes less accurate. 

For example, a crude estimation of the average energy loss per col- 
lision yields values of 25.9 eV and 2.9 eV' for ionization and excitation 
collisions respectively for an electron at 32 eV, giving a AE/E of 
approximately unity, far too large for continuous slowing. In contrast, 
a similar estimation can be made for Coulombic collisions employing the 
expression for the electron collision frequency, 

~ ” C^1Y) 

jf'P 

where is replaced by the velocity of an electron at 32 eV and the 

Coulomb logarithm is approximately 2.8. This yields a collision frequency 
of M.S X lO^^sec or '\'0.044 eV lost/collision. Then AE/E is ''^0.001 
which is quite consistent with the continuous slowing assumption. Thus, 
as seen from Table 1, the results are expected to be more accurate for 
high temperatures where Coulombic collisions contribute a large fraction 
of the energy loss rate. 
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TOTAL 

COULOMBIC 
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0 


u" 


COLLISIONS 

IONIZATION 

COLLISIONS 

EXCITATION 

COLLISIONS 

IONIZATION 

COLLISIONS 

EXCITATION 

COLLISIONS 

eV 


2.26 

.14 

.60 

1.45 


.01 

.06 

826.2 


4.23 

.24 

1.13 

2.70 


.03 

.13 

179.3 

5000®K 

5.53 

.36 

1.28 . 

3.68 


.03 

.18 

57.0 


5.86 

.45 

1.14 

4.04 


.02 

.21 . 

26.0 


.60 

.60 


. ^ - 


- 

- 

3.0 


2.13 

.41 

.43 

1.03 


.05 

.21 

826.2 


3.94 

.72 

.80 

1.92 


.10 

.40 

179.3 

6000‘*K 

5.24 

l;04 

.91 

2.62 


.09 

.58 

57.0 


5.71 

1.29 

.81 

2.87 


.06 

-.68 

26.6 


1.40 

1.40 

- 

- 




3.0 


2.17 

.81 

.25 

.60 


.09 

.42 

826.2 


4.00 

1.40 

.47 

1.12 


.19 

.82 

179.3 

7000°K 

5.42 

2.00 

.53 

1.53 


.19 

1.17 

57.0 


6.11 

. • 2.45 

.47 

1.68 


.13 

1.38 

26.6 


2.17 

2.17 

- 

- 


- 

- 

3.0 


2.23 

1.13 

.11 

.27 


.13 

.59 

826.2 


4.06 

1.93 

.21 

.50 


.27 

•• 1.15 

179.3 

8000“K 

5.57 

■ 2.75 

.24 

.68 


.26’ 

1.64 

57.0 


6.42 . 

3.34 

.21 

.75 


.18 

i.94 

26.6 


2.52 

2.52 


- 



- 

3.0 


Table 1. Energy loss rates listed as a function of energy for various temperatures. 
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E . Parametric Results 

1 . Temperature Variation 

The results o£ calculating the distribution function according to 
the analytic prescription of Eq. (10) for various, temperatures are plot- 
ted in Fig. 5. The most noticeable effect of temperature variation is the 
increased magnitude of the deviation of the high-energy tail from the 

Maxwellian with decreasing temperature. This effect is directly traceable 

f331 

to the energy dependence of the average fission cross section employed. 
As the .plasma temperature and the corresponding average energy of the 
thermalized neutrons is decreased, the fission cross’ section increases, 
■ultimately, yielding a larger nascent electron source rate (see Fig. 6) . 

The degree of deviation from a Maxwellian as well as the energy range for 
which the non-Maxwell ian behavior is dominant is thereby enhanced with 
decreasing temperatures. The point of intersection of the high-energy 
tail with the Maxwellian distribution denotes the lowest energy for which 
Eq'. (10) is valid. 

The effect of temperature variation upon the slopes of the. high- 
energy tail is extremely subdued over most of its energy range.. Only at 
the lowest energies, i.e., at the intersection of the high-energy tail with 
the Maxwellian, is there any noticeable difference. Examination of the 
energy "loss rates in Table 1 reveals a partial explanation for the be- 
havior of the slopes. The energy loss rates for the range of tempera- 
tures considered are more disparate at low energies. The inelastic cross 
sections fall off drastically at low energies, accounting for the low 
energy behavior of the energy loss rate while the temperature dependence 
of the density results in the energy loss rate being -nearly independent 
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of temperature, at high energies. 

A second factor influencing the slopes is the secondary electron 
source, - Not. only are the magnitudes of the sources different for the 
various temperatures considered but the shapes of the source distributions ‘ 
are different at' low energies. The latter is caused by variations in the 
fraction of neutrals present in the plasma compound by difference in 
electronic stracture between neutral uranium and singly ionized uranium. 
(This effect is observed only at energies near the threshold for ioniza- 
tion.) 

2. Flux Variation 

The effect of neutron flux variation upon the distribution function 
is considerably less complex than the effect of temperature variation. 

Under the plasma conditions studied, the bulk of the thermalized electrons 
is the result of the high plasma temperature. The neutron flux does not 
alter the thermalized densities by more than 1/10% from -the normal Saha 
values. Then, the only effect a change in the neutron flux level can 
produce is- a change in the production rate of high energy electrons. 
Therefore, according to Eq. (10), the high-energy tail is directly pro- 
portional to the neutron flux level j which is consistent with the results 
in Fig. 7 where the high-energy tail calculated at one flux level' is 
simply a scaled vertical translation of the tail at a different flux 
level . 

3. Cross-Sectional Dependence 

A vital aspect of the interpretation- of data is its credibility. The 
largest inaccuracy existing in the calculation of the distribution function 
lies with the uncertainty in cross-sections. Since no ej{periraental data 
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Fig. 7. The effect of varying the neutron flux upon the distribution 
function for a constant temperature of SOOO^K. 
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exists for uraniim, a hydrogenic model formed the basis for calculating 
the necessary cross-sections (see Appendix B) . A comparison of the 
cross sections to measured values for cesium and helium reveals the 
uranium cross sections to fall between them, but cl9ser to helium. 

However, due to the similarity in electronic structure of cesium and 
uranium, the uranium cross sections would be expected to lie closer to 
cesium than helium. To investigate this, another set of cross sections 
are obtained by doubling the inelastic cross sections. If the doubled 
cross-section set is used in the calculation of the distribution function 
in conjunction with the Coulombic energy loss rate predicted by the mii- 
fied theory (see Appendix A), the possible errors generated by inac- 
curate slowing theory can be gauged. 

The result of just such a calculation is compared with a calculation 

with the unadjusted cross-section set and the Fokker-Planck slowing theory 

in Fig. 8. Fortunately, the differences indicated are not large, e.g., a 

maximum deviation of 4% is observed at '^'20 eV for 5000®K. Insight as to 

the reason for the differences can be gained from the energy loss rates 

appearing in Table 2 which were employed in this calculation and those 

in Table 1 for the previous calculations. Due to the doubling of the 

inelastic cross sections, both the source and the inelastic energy loss 

dE 

rates are doubled. Since f(E) “ factor of two is cancelled, 

provided the inelastic events dominate, as they do at 5000“K. At 8000®K, 

I 

however, the increase in the source rate, a result of the ionization cross 
section being doubled, is not totally compensated by an increase in the 
inelastic energy losses. This occurs because of the- relatively large 
contribution by elastic collisions at higher temperatures. 
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Fig. 8. The distribution function (-) versus energy compared to repeated 
calculation of the distribution with the cross-section set 
doubled (-•-). 
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eV 


4.40 

.14 

1.20 

2.90 

.. 02 

.13 

826.2 


8.21 

.24 

2.26 

5.40 

.06 

.25 

179.3 

5000'*K 

10.69 

.36 

2.55 

7.36 

.06 

.36 

57.0 


11.29 

.47 

2.28 

8.07 

.04 

.42 

26.6 


.87 

.87 

- 

- 

- 

- 

3.0 


3.33 

1.13 

.22 

.54 

.26 

1.18 

826.2 


6.24 

1.97 

.42 

1.00 

.55 

2.31 

179.3 

8000° K 

8.54 

2.88 

.47 

1.37 

.53 

3.29 

57.0 


9.80 

3.64 

.42 

1.50 

.37 

3.88 

26.6 


5.75 

5.75 

— 


— 

- 

3.0 


Table 2. Energy loss rates listed as a fimction of energy for various temperatures 
with the inelastic cross sections doubled and a unified theory treatment 
of the Coulombic collisions. 


w 











32 


The errors projected in this section are indicative o£ those antici- 
pated to appear in the calculation of the distribution function. The 
inelastic cross sections display a shape characteristic of other elements, 
most notably cesium, but are low in magnitude over all energies of in- 
terest by an estimated factor of two. The unified theory expression for 
the Coulombic energy loss rate more accurately depicts collective and 
binary interactions than the corresponding Fokker- Planck expressions, 
yet the difference is not so large as to discredit the Fokker-Planck re- 
sult. From these results, it is seen that the anticipated inaccuracies 
in the inelastic cross sections and the Coulombic energy loss rate are 
compensating inaccuracies, i.e.,. the inaccuracy of the distribution is' 
less than the inaccuracies associated with either the inelastic cross- 
section set or the Coulombic energy loss rate. 



33 


CHAPTER II r 
MONTE CARLO TECHNIQUE 


A. Introduction 

In the last chapter, it was pointed out that the assumption of con- 
tinuous slowing is questionable, particularly at low energies. To treat 
the problem in a more precise manner, an improved treatment of inelastic 
collisions is required. Due to the increased complexity of the present 
treatment of the slowing process compared to the continuous slowing 
treatment, a Monte Carlo simulation was selected and is described here. 

A straightforward approach is to follow the electrons via an analytic 
prescription for Coulombic collisions' for a time equal to the inverse of 
the inelastic collision frequency. At that time they suffer a discrete 
inelastic event which is treated by normal Monte Carlo techniques and 
then the process is repeated. The present method is an improvement over 
this technique in that provisions are made for the variation of the in- 
elastic cross section between such events. Of course, the energy at 
which the inelastic collisions occur as well as the energy loss suffered 
are chosen in a random fashion according to the appropriate probability 
distribution. 

Bj Derivation of Governing Equation 

The results of the previous section indicate the need to relax the 
assumption of continuous slowing down for inelastic collisions. An ap- 
propriate equation may be derived from Eq. (1), namely: 


SceME + o, = a- 

Jin Jo 


( 18 ) 
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E E + dE 



Fig. 9. Illustration of collisional processes used in the development of 
particle balance in energy space. The processes are: 1) inelas- 
tic scattering out of the interval dE about E 2)' elastic CCoul- 
ombic) scattering out of the intei*val 3) secondary electron 
production in interval 4) elastic scattering into interval 
5} inelastic scattering into interval . 
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where, as before, recombination is neglected. The slowing terms in 
energy space are readily obtainable from the diagram in- Fig. 9. 

The slowing down processes labeled 2 and 4 appearing in Fig. 9 are 
characteristic of the continuous slowing down model previously employed, 
but here it will be restricted to the elastic scattering component. 
Hence, 


$(£■) A£ +f(E + dlE) 4x1 +■ 'ih - ‘f-'f 

where the primed terms indicate inelastic scattering. Expressions for 

and may be obtained with the Gryzinski^^^^ type energy' transfer 

differential cross section ^ (E,e) introduced in Eq. (13). The terra 
q'out (19) represents scattering via an inelastic process (repre- 
sented by arrow 1 in Fig, 9) from the energy interval dE about E (the 

shaded region (B) in Fig, 9) into any energy interval dE below energy E 

(region A in Fig. 9). Mathematically, this can be written as. 




(K/ti. 


= n 





l <i& 

die 


(20) 


where e is the energy lost per inelastic scattering event, n is the 
density of the background species able to participate in the particular 
event under consideration, and the relative velocity is approximated by 
the velocity of an electron of energy E, i.e., . The range imposed 

’ wig 

upon the energy lost e varies from 0 to E which is determined by the 

max 

process involved • Since there are a number of ionization and excitation 
processes competing in the slowing process, a sum over these processes is 
necessary, i.e.. 
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r / — r ^ 1 

Similarly, the scattering process labeled 5 represents electrons of an 
energy greater than E (i.e., electrons from region C) which- scatter into 
the energy interval dE about E, i.e.. 



Substitution of Eqs. (21.) and (22) into Eq. (19) yields Eq. (23): 

SCE)<t£: +f(.E:+<iE)4E) 

c£t /e+<te 






y. {nr j ~y FCE+£ ) de'iE+tX) 


"le ct£; 



i dE 


■ J ^ -X 

= t [n.j^ yt <££ 


+* 


-PcE-1 




£ 


(23) 


The source term S(E) appearing in Eq. (20) may be decomposed into a 
series as was done in Chapter II. The first terra S^ in the series repre- 
sents the nascent electron source while additional higher order terms 
depict the various generations comprising the avalanche of secondary elec- 
trons. A concise expression for the total secondary electron soxirce may 
be derived from considerations of the ionization process labeled 3 in 
Fig. 9. The secondary electrons born in the energy interval dE about E 


are the result of ionization collisions in which an incident electron 
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loses an energy E + where E is the secondary's kinetic energy and 
is the energy necessary for the secondary electron to overcome the ioniza- 
tion potential of the target uraniiira atom. Due to the indistinguishability 
of electrons,, the least energetic of the resulting pair is defined as the 
secondary electron and the other, more energetic electron is defined as 
that electron which was termed the incident electron before the collision 
occurred. Then, integration of the ionization reaction- rate over the 
energy range of electrons capable of generating a secondary electron of 
energy E followed by a summation over the various bound electronic states 
which can participate in secondary electron production yields the follow- 
ing expression for the total electron source rate: 




'• *■ ^ '■ 

The calculation of the nascent source term S^' is outlined in Chapter II. 


Then, 


S„(E)d£ 
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The governing equation .can be obtained from Eq. C25) by invoking the. 
fundamental theorem of calculus, i.e.. 


- fe 


ieT'-ZHilS}. A(9'CE+£,e)fc£:+t> dLs- 

t • M) ' >^e 
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*- CE+Vi) ^ *^e cJLf-c 
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x^2.E 


mOL.X; 


X-niJ ir f^£ia-P(E) 


(26) 


The inelastic collision terms in Eq. (26) have been derived elsewhere. 
However, their appearance with the continuous slowing down treatment is 
unique . 

C. Simulation 

1. Distinction Between Collisional Processes and Their Treatments 
a. Coulombic Collisions 

Equation' (26). is not .amenable to an analytic solution as was Eq. (7), 

but it does lend itself- to a novel method of solution involving a Monte 

Carlo simulation that integrates both analytic and random walk techniques. 

An analytic description of the Coulombic collisions is employed in the 

tracing of the histories of electrons. An individual electron is permitted 

to evolve for a time At as prescribed by the following equation. 

,t+ At 


5 


dLt = AE 

iCou\omblc 
• coir.s^on 


(27) 


where ^ 
at 


Coulombic ^®P^®sents the electron energy loss rate due to Coulombic 
collision 



39 


collisions and AE is the energy lost during the time At. Then, if 
an electron has an energy E associated with it at time t, at time t + At 
its energy is E - AE. The inelastic collisions are then superimposed on 
the Couloinbic slowing down in a discrete fashion as described below, 
b. Inelastic Collisions 

i. Choice of Collision Energy 

Since the inelastic collisions are less frequent than the Coulombic 
collisions, they are superimposed in a random walk fashion. The distance 
of the walk is prescribed by a probability distribution dependent upon 
the test particle's energy, the inelastic cross-sections, and the 
Coulombic energy loss rate. Such a function can be obtained by first 
examining the probability PCx) of a collision occurring in an infinitesimal 
distance dx about x measured along the flight path of an electron. This 
probability is a product of the macroscopic cross section Z for inelastic 
collisions and the length of the interval dx, i.e., Zfx) dx. The 
functional dependence upon x is included as a reminder that the inelastic 
cross section depends on energy which in turn is dependent upon the dis- 
tance of travel within the slowing medium. The density of the target 
particles is assigned to be constant within the medium. 

The probability of traveling a distance x without a collision is the 
ratio of the intensity of a beam of test particles displaced a distance x 
to the initial intensity at x = 0, i.e., lCx)/I(0). The attenuation of 
such a beam is governed by the following equation: 


— z, CTO I C7-) AtL 


(28) 
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Integration of Eq. (28) yields, 

Xc%) 

Tio) = \ZaUt' ) 

Then, the probability of the first inelastic collision occurring in dx 
about X is 

Since the problem is to be solved in energy space, an equivalent probability 
of the first inelastic collision occurring in the energy interval dE about 
E is desired, or 




C31) 


thereby implying that the Jacobian necessary for such a transformation would 
dE 

be ^ which can be related to the energy loss rate by. 


Coelom bit 

d% u 


C32) 


Then, the transformation of Eq. C30) into the energy variable E yields 
Eq. (33): 


FpCE,)dE = LcE') 


2£ 
1 v«e 


Jl c* I ^ 

^ j coU;sl<in5 ^ 

it 


■liM, 


irSc^E (33) 

rt f~ 

CL^ j<!o>Us? 0»lS 

where P^CE) is the probability that the first inelastic collision will 
occur at energy E if the electron is injected at energy E^. 
ii. Determination of Energy Lost 
Once it has been established that a collision occurs at energy E, 
the amount of energy lost must be determined. This too is prescribed by 
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a probability distribution function. The probability of a particle which 
collides at energy E losing energy e may be obtained by normalizing, the 

Arf 

differential energy transfer cross section (E,e) to the microscopic 
cross section ffCE) for that process, i.e., 

Ae C34) 

AL 6'Ce) 

The determination of the process, i.e., which of the uranium species, 
collision types, and atomic levels are involved, must preceed the energy 
'loss calculation (see Appendix C for the algorithms employed) . 

2. Merging, of Treatments 

a, No Inelastic Collision 

The manner in which all of the above aspects of collisions are in- 
corporated to yield a complete description of the slowing down process 
will now be illustrated with an example. Let us begin by assuming an 
electron is born at energy which corresponds to a time t^ on the energy 
vs. time plot on the right hand side of Fig. 10. This plot, generated 
according to Eq. (27), represents an electron's energy as a function of 
time as it slows doim solely due to Coulombic collisions in an inter- 
mediate energy range well above E^. Since the introduction of source 
particles must occur frequently enough to approximate continuous inter- 
jection, an individual electron is only permitted to evolve for some small 
time interval At. The result of one such period is to let the electron 
follow the slowing down curve to the point corresponding to time t^^ = t^ + At 
and energy . A second period leaves the electron at time ^2 ~ + 2At 

and energy This process continues until that period in which the 
electron falls below the energy E^ for which Eqs. (7) and (13) become in- 



Fig. 10. Scheme for electron simulation. Figure contains plots of 

energy versus time as predicted by Eq. C27) , of the collision 
prqbability P , versus collision energy from Eq, C37) and of 
the probability of losing energy AE in an inelastic collision 
Pac versus final electron energy from Eq. (38). 
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valid, or an inelastic collision occurs, 
b. Inelastic Collision 

The decision for the existence of an inelastic collision must be 


made at the beginning of the electron’s history, i.e., at E^. It is based 
upon the probability of an electron traversing the energy range to the 
lowest valid energy considered E.p, or the non-collision probability: 


' rr rtJfc. UotUsIonS 


At k‘ 

This probability is by necessity equal to one minus the integral of the 


probability of a collision occurring at energy E, Eq. (33) , with the limits 


of integration being from E^ to E^. Comparison of with a random num- 
ber chosen from the range 0 to 1 completes the decision process. If the 


random number is greater than Pj^^, then an inelastic collision must occur. 
Conversely, if the random number is less than Pj^^, the electron will not 
suffer an inelastic collision. 


For the sake of thoroughness, assume the electron collides. Then, 

the energy at which the collision will occur must be selected. Normally, 

the energy would be randomly selected from the inverse distribution of the 

probability of colliding at energy E. However, if the probability dis-- 

tribution is too complex to invert, as it is here, a form of the rejection 
(34) 

technique must be employed. This algorithm begins by mapping the 
probability distribution onto a rectangle of unit area. The prescription 
for obtaining an acceptable candidate for the collision energy E involves 
the selection of two random numbers, r^ and X 2 > where r^^ represents an 
evaluation of the probability distribution and the other, r 2 > is a candi- 
date for E. If r^ is less than the distribution evaluated at then T 2 
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is accepted as the collision energy. However, if the converse is true, 
the candidate is rejected and the process is repeated until an acceptable 
candidate is found. CThis process is analogous to throwing darts at a 
rectangle with the coordinates of the impact point corresponding to r^^ 
and r 2 « When the impact point falls below. the probability distribution, 
the corresponding energy is taken to be the collision energy.) 

As applied to the calculation of the energy of the next collision, 
the technique commences with the following prescription for the renormal- 
ized collision probability: 


P - 




(36) 


where the expression for is obtained from Eq. (33) and the energies 
^col ^o respectively, a candidate for collision energy and the 

initial electron energy. Inherent to the success of the algorithm is 
that the raaximvun value of the probability distribution is readily ob- 
tainable, i.e., the maximum must occur at E . Then, 
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A plot of P^q 2 "''^®^sus appears in Fig. 10. 

The choice of collision energy is completed by first randomly choosing 
a candidate from the energy range E^ to E^. Next the corresponding 
= R is calculated and compared with a random number in the interval 
0 to 1. If the random number is less than R, the candidate energy E^ is 
accepted as the site of the collision. If the random number is greater 
than R, the process is repeated until an acceptable candidate is obtained.. 
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Continuing the present example, the collision energy is assumed to 
be E_. Then the electron slows from energy E to E_ via Coulombic col- 

O 0 o , 

lisions. In the third time period, the period in which the -collision 
occurs, the electron would have proceeded from time t2 and energy E2 to 
time t^ = t2 + At and energy E^ if a collision had not occurred. The 
occurrence of a collision does not alter the fact that the electron must 
interact via Coulombic collisions for the full time period At (the in- 
elastic interaction time is negligible compared to At) . Then, after the 
occurrence of the inelastic collision, Coulombic interactions are to be 
resumed for a time At’ in order to complete the evolution of the electron 
for the period At. 

The detennination of which of the background species is involved 
and also which of the possible ionization and excitation events will be 
involved in the collision must precede the calculation of the energy' lost 
as a result of the collision. The species selection is accomplished by a 
comparison of a random number (henceforth in the discussion all random 
numbers are assumed to be evenly distributed from 0 to 1 ) to the ratio 
of the macroscopic cross section of a species to the total macroscopic 
cross section in the usual Monte Carlo fashion. Similarly, the type of 
collision is chosen by comparison of a random number with the ratio of 
the microscopic cross section for a species process to the total 
microscopic cross-section. 

For economical reasons, the energy lost in an excitation collision 
is approximated as the excitation energy. This approximation is a 
reasonable one because the possible energy losses range from the excita- 
tion energy of the process considered to the excitation energy of the net 
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,state affording a larger energy transfer. Furthermore, the difference in 
energy of two adjacent excitation levels is usually small compared to the 
excitation energy (see Table 4) and the probability of an energy transfer 
event increases as the amount of energy transferred decreases-. 

In the case of an ionization collision, the amount of energy lost 
by the electron at energy is also determined by a rejection technique. 
In this case, the renormalized probability P^ p is 



(38) 


where P^(AE) is the probability of losing AE energy through an ionization 
collision and is obtained from Eq. (34) . The energy Ej^ represents the 
minimum energy which can be lost and is non-zero because of the presence 
of a threshold energy necessary to initiate the ionization process. 

In Fig. 10, a plot of P^ is presented as a function of the final 
energy (rather than the energy transferred) . As before, a candidate 
is chosen as the .energy of the electron after the collision and the cor- 
responding renormalized probability P^(Eg) is subjected to the acceptability 
criterion. Assuming R' to be larger than the random number chosen, the 

tracing of the electron can proceed from energy Ej.. 

o 

At this point a decision is made whether or not the electron will 
incur further collisions and where in energy the next collision will occur. 
This is done as previously described. After performing these tasks, an 
attempt is made to permit the electron to proceed a time At* such that 
the period will be completed. The electron is advanced to energy E^ 
corresponding to time t^ = tp + At*, provided a collision does not occur at 
an energy greater than E^. Should another collision occur in the same 
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period, the previous process is repeated until the electron has been per- 
mitted to evolve for the duration- of the time period At, and so the 
evolution continues for additional periods, moving in steps .of At along 
the Coulorabic slowing down curve, to the next collision and eventually 
past the threshold energy of validity E^. 

3. Computational Time Reduction 

Five additional techniques are utilized to provide a substantial re- 
duction in computational time without altering the accuracy of the code. 
The techniques are described below in the order of their effectiveness. 

The first technique is a unique, new method called convergence 
propagation. A necessary cri-terion for its applicability is that f(E) 
is dependent upon f(E’ > E) (rather than f(E’ < E)). Then, capitalizing 
upon this condition, the method is the dynamic expansion downward in 
energy of the energy region over which the distribution is simulated. 

The expansion of the simulation region occurs only when convergence has 
been obtained in the current region. The savings in computational time 
is realized by not having to simulate the distribution below energy E 
until the wave of convergence has arrived at E (for a more detailed dis- 
cussion of this technique, see Appendix C) . 

A second technique involves the fitting of frequently' evaluated, 
complex functions with a number of quadratic equations, each valid in a 
unique subinterval of the dependent variable's range. The coefficients 
of these equations are determined by a cubic spline algorithm (from the 
IMSL subroutine library). The savings are substantial, and better than 
single precision accuracy is easily obtained for the functions. 
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A third technique is strongly coupled to the convergence propagation 
technique in that an optimum time step is chosen for each energy group in 
the convergence scheme. The time for ah electron to traverse the yarious 
energy groups varies by four orders of magnitude in the present problem. 
Thus, it becomes essential to. gear the frequency at which convergence is 
checked, i.e., At, to the group transit time of the current converging 
group as the convergence wave propagates towards lower energies. Devia- 
tions in the specification of At as the transit time are permitted by 
examining the ratio of the nvimber of electrons in the group to the number 
of source electrons introduced. If the source is the major input into 
the group, then At is shortened from the transit time. Conversely, if 
the source is not the dominant input, then a larger time is used to per- 
mit sufficient collisions to occur in order to get better statistics on 
the- input into the group due to inelastic collisions. By this scheme, 
the convergence check will be made as soon as a significant change has 
been made in the -distribution function. 

A fourth technique relies upon stacking source particles in energy 
and staggering their associated time of introduction. At the time of 
particle introduction or replenishment, m particles, where m = k*£,(m,k,il 
are integers), are introduced into an energy interval, where 5. particles 
are assigned the same injection energy, there being k such energies. To 
avoid the m particles behaving as if only k particles had been introduced, 
each is assigned a unique collision energy and associated time of injec- 
tion. Since particles are replenished with a periodicity of At, an 
assigned injection time t must lie randomly within the interval 
T^ - At < t < T^ is the time at which the current replenishment occurs 
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according to the internal clock. Many facets of the time consuming cal- 
culation of the collision energy are thereby retainable for the other 
A-1 particles at the same energy. 

A fifth technique involves the utilization of an initial guess. The 

proximity of the guess to the ultimate answer determines the efficiency 

of this technique. However, the effectiveness of this technique as re- 
f23") 

ported by Wang^ ^ was not realized in the present Monte Carlo code. This 
is attributed to the effectiveness of the convergence propagation tech- 
nique. 

D. Test Run 

Immediate questions arising after the development of a Monte Carlo 
code concern the accuracy and the rate of convergence to the solution of 
the proposed problem. The criterion for convergence involves various 
conditions such as: consistency, approximity to an experimentally ob- 

served quantity or independent calculation, and stability of the solution. 

Each of these tests are stringently applied to an arbitrarily chosen case 
in the ensuing paragraphs. 

16 2 

The case in question is for a neutron flux of 2 x 10 neutrons/ (cm"^ -Sec) 

and a plasma temperature of 5000"K. A sample of the results is shown in 
Fig. 11. The figure is a composite of four graphs, each depicting the 
high-energy tail at successive times as recorded by an internal clock. 

Also appearing in each graph is the tail of a Maxwellian distribution cor- 
responding to the electron density and temperature, and the corresponding 
analytic solution for the high-energy tail. As in Chapter II, the total 
distribution function in either the Monte Carlo or analytic calculation 
is the union of the Maxwellian and the corresponding high-energy tail. 
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Fig. 11. Snapshots of Monte Carlo distribution at four, different 

times. Also displayed are a Maxwellian distribution (light 
line) and the analytic solution (dark line) , 
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The relationship between the analytic and Monte Carlo result is atypical 
of other temperatures and shall be explored in detail in Section F of 
this chapter. . For the present, the analytic result will serve solely as 
a .standard, for comparison. The final result to be reported later is an 
average of a dozen such snapshots of the distribution function. 

1 . Consistency 

Since a computed result may seem plausible and yet be erroneous, 
especially with large computer programs such as the present Monte Carlo 
code, the code must be established to be free of programming and logic 
errors. This is accomplished' by comparing, for consistency, computer 
calculations with independent hand calculation. 

The aforementioned computer results are not those of the distribution 
fiinction, but rather additional, supplemental results relative to the dis- 
tribution function, recorded at the time of each snapshot. A sample of 
these results, corresponding to the last snapshot of Fig.- 11, appears in 
Fig. 12. Each 'entry in Fig. 12 is briefly described below. 

The energy lost by all the electrons within an energy group centered 
at the tabulated energies is recorded, providing information regarding the 
slowing- mechanisms. The presence of zeros in the ionization and excita- 
tion energy loss rate columns indicates there were no collisions of this 
type during the last time period At. 

Also, the W-value (the energy lost per ion pair formed) is calculated 
during the same time period. This particular value is larger than the 
■average value based upon results at eleven other time periods, i.e. , 

160 eV per ionization event. The. importance of this particular result 
cannot be fully realized as there are no experimental measurements with 
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which to compare and the result is plagued fay a large statistical un- 
certainty. However, the result is consistent with the expectation that 
xt be larger than the ionization potential of 6.2 eV for neutral uranium. 

The time elapsed on the internal clock is presented in comparison 
with an estimate of the time necessary to achieve a converged solution 
by a straightforward Monte Carlo simulation. If the elapsed time is 
larger than the estimated time, then the electron slowing calculation is 
grossly in error even if the initial guess is very poor. Then, a value of 
3.6 for the ratio of the estimated time to the elapsed time is reassuring 
concerning the validity of the electron slowing. This factor of 3.6 is 
also a measure of the efficiency of the convergence propagation technique. 

The data necessary to determine if particle balances exist for the 
lowest energy group and the entire ensemble of energy groups appears in 
Fig. 12. A particle balance to within a tolerance of ±15% is imposed upon 
the current lowest group as a condition for advancing to the next group 
in the convergence propagation scheme. If the lowest energy group complies 
with this condition, then a snapshot of the high-energy tail is recorded. 

No provisions are made to guarantee a global particle balance. Then, the 
observance of such a balance to within ±1% insures the validity of the 

propagation of convergence technique and the accuracy of the integral of 
the distributioji function* 

Further consistency checks can be performed with the data in Fig. 12 
For example, the total Coulombic energy loss rate at 25.1 eV divided by 
the number of particles within the group, f CE) dE (where it is correctly 
assumed that all the electrons in the group suffer this collision-type) 
yields an energy loss rate of .41 ergs/sec. The average value obtained 
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for a dozen snapshots is .45 ergs/sec. These values are within ten per- 
cent of the calculated value, in Table 1, of .448 ergs/sec. 

For inelastic collisions, the energy loss rate is tabulatsd as the 
result of two- separate algorithms; one establishes the collision frequency 
while the other determines the energy loss. A separate run specifically 
designed to check the average energy loss produced results very similar to 
those calculated in Section D of Chapter’ll. The mean value of the energy 
loss for excitation collisions, 2.98 eV, obtained from the code is com- 
-parable to an average value of 2.9 eV calculated in- Chapter II. For 
ionization events, the code generated values of 23.4 eV compared to an 
average value of 24.5 eV also calculated' in Chapter II. 

The accuracy of the collisioiji frequency algorithm can also be 

I 

verified by comparing the niiraber of excitation collisions generated by the 

code with the number of ejq)ected excitation collisions as calculated below. 

Dividing the excitation energy loss rate at 25.1 eV of Fig. 12 by an 

average energy loss of 2.98 eV yields a value of 3.43 x 10^^ excitation 

collisions per second. Alternately, since there are 8.3 x’ 10^ electrons 

at 25.1 eV and the collision frequency is 8.7 x 10^ coliisions/sec, 

IS 

7.2 X 10 excitation collisions are expected per second. These results 
are in reasonable agreement, and. the discrepancy can be attributed to 
statistical fluctuations as the average excitation energy loss obtained 
from the Monte Carlo code is 3.2 jc 10^ ergs/sec compared to the value of 

4 

1.63 X 10 ergs/sec reported for this particular snapshot of the dis- 
tribution function. From these checks, the code' can be concluded to be 
a consistent and accurate representation of the physical processes- in- 


volved. 
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2. Approximity to Independent Result 

Confidence in the Monte Carlo calculations can also be gained by a 
comparison with the results -o£ the earlier analytic model. The two re- 
sults are in exceedingly good agreement in Fig. 11. However, such excep- 
tional agreement is atypical (see Fig. 18). Nevertheless, the assumption 
of continuous slowing down, which is the basis of the analytic model, is 
not too unreasonable under the plasma conditions considered here. There- 
fore, the approximity of the two models* results is expected and its ob- 
servance (see Fig. 18) reaffirms our confidence in the Monte Carlo calcu- 
lations. Furthermore, the existence of a display of common trends (see 
Section E) and a predictable disparity (see Section F) reinforce the 
acceptability of the Monte Carlo solution. 

3. Stability 

The final criterion for the acceptability of the solution is its 
exhibition of stability. The stability of the "converged’* solution is dis- 
played in Fig. 11 over a short period of time while a longer term time 
history of the distribution evaluated at a single energy’ is shown' in Fig. 
13. In both instances, random oscillations about an average value- are 
observed. The display of this type of behavior in conjunction with the 
observance of a global particle balance suggests that the solution has 
converged. However, these criteria for convergence .are not acceptable 
until they have been demonstrated to correctly predict the solution to be 
stable for extended periods of time (at least several multiples of the 
Coulombic slowing down time). 

The necessity of observing the distribution then for tens of nano- 
seconds imposes a severe financial strain. Thus, a less costly, but 



Fig. 13. The history in time of an individual point fCE) on the Monte 
Carlo curves of Fig, 11 (A's) 


30. 7r 



A 


STANDARD DEVIATION = 2.85x10^ electrons/(cm^-eV) 
AVERAGE VALUE =2.49 xio"^ electrons/(cm^-eV) 

DATA FOR ENERGY GROUP 1 1.6 eV 


Ln 

ON 



57 


equivalent scheme of sustained observation, was ultimately employed. The 
scheme relies on the knowledge that a necessary condition for terminating 
observations is the elimination of any correlation in time of the last 
observation with the initial observation. In Monte Carlo calculations, 
this is equivalent to using a random number string of infinite period 
or several strings of large but finite period. Then, by repeating the cal-^ 
culation several times, each time using a different random number string, 
sufficient data will exist to determine if the solution has relaxed into 
a stable configuration. 

Such tests were performed for a distribution at 8000“K. This partic- 
ular temperature was chosen because the final distribution differs most 
from the initial guess, the analytic solution. The results of a dozen 
observations taken in six separate runs of the code, with different ran- 
dom number strings for each run, are swnmarized in Table 3. Each average 
of the dozen observations for the six runs falls within a standard devia- 
tion of the average of all seventy-two observations, except for the values 
at the two lowest energies. The reason for the bad statistics at low 
energies lies in the choice of the machine particle distribution Csee 
Appendix C) . Since these two distribution points are past the intersection 
of the Maxwellian and are not meaningful, their behavior is irrelevant. 

The observance of oscillations in the standard deviation are the result of 
the coarseness imposed upon the calculation through the number of simula- 
tion particles utilized, combined with the fact that an electron will, on 
the average, suffer three inelastic events while slowing from 1 keV to 
3.0 eV for an 8000"K plasma. 
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These tests were repeated at 5000®K, but with only three separate 
runs (and three random number strings) . As before, an average of the 
dozen observations recorded in a single run falls within a standard 
deviation' of the average of all thirty-six observations. Furthermore, 
at 8000®K and SOOO^K, the standard deviation calculated in each run 
approximately equals the standard deviation of all observations at the 
same ten 5 )erature. Then, a single run yields sufficient data from which 
to conclusively generate the converged solution. The converged solution 
is the average of the dozen observations recorded in a single run to 
within a standard deviation, also calculated in the run (see Table 3 for 
typical values of the standard deviation). 

E. Parametric Studies 

1 . Temperature Dependence 

The resulting distribution functions of Monte Carlo calculations 

are displayed in Fig. 14 for various temperatures at a constant neutron 
14 2 

flux of 2 X 10 neutrons/ (cm -sec) . The overall trends indicated by the 
Monte Carlo results are quite similar to the earlier analytic solutions- 
in Fig. 5. The distribution function appears to be dependent upon tem- 
perature variations predominantly through the normalization of the high 
energy tail, except for a slight change in slope at low energies. 

2. Neutron Flux Dependence 

Similarly, the distribution functions calculated by the Monte Carlo 
code reflect the same trends in parametric variations of the neutron flux 
as do the analytic results, namely, the distribution function is directly 
proportional to the neutron flux level. This is amply illustrated in 
Fig. 15 for a plasma temperature of 5000“K. 



59 


DISTRIBUTION FUNCTION, electrons/ (cm^-eV) 
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Table 3. Statistical data on 'Monte Carlo results for 8000“K and a 
neutron flux of 2 x 10l4 neutrons/ Ccm^-sec) based on a 
dozen observations each of six runs with different random 
■ number strings . 
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Fig. 14. Temperature 'dependence of distribution function is exhibited 
for a neutron flux of 2 x 10^4 neutrons/ (cm^-sec) at the 
temperatures of 8000*K, 7000“K, 6000®K and 5000“K. 
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3. Comparison o£ Coulorabic Energy Loss Rate Models 
Previous to this, all distribution functions have been calculated 
with a traditional Fokker-Planck expression for the energy loss rate for 
Coulombic collisions. Due to the small number of electrons in a Debye 
sphere, collective interactions cannot be overlooked.. The distribution 
function has been recalculated with a Coulombic energy loss rate from the 
'unified .theory which incorporates both binary and collective inter- 
actions into a single theory. The results are presented in Fig. 16 for 

14 2 

a neutron flux of 2 x 10 neutrons/ (cm -sec) and temperatures ranging 
from 8000®K to 5000'’K. Also, for completeness, the range of electron 
energy has been extended closer than previous results to the maximum 
energy at which a nascent electron can be bom (-2.1 keV) 

Again, the general trends of the previous calculations are still 
preserved under a change of expressions for the energy loss rate. How- 
ever, the absolute magnitude of the high-energy tail is affected by. the 
change. This is illustrated in Fig. 17 .for a neutron flux of 2 x 10^“^ 
neutrons/ (cm -sec) and a temperature of 8000°K. (This represents the 
"worst” case since at higher temperatures, the Coulombic energy loss rate 
.comprises a larger fraction of the total energy loss rate than at any other 
temperature considered.) In general, the energy loss rate for the unified 
theory is approximately 1.5 times larger than the Fokker-Planck energy 
loss rate, due to the additional slowing mechanisms (collective inter-r 
actions and hard impact collisions) considered in the former energy loss 
rate. This factor decreases the distribution function by approximately 
a factor of 2/3 as predicted by Eq. (10) . 
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ENERGY, eV 

Fig. 16. Temperature dependence of distribution function is repeated 
as in Fig. 14, but this time using the unified slowing 
treatment of Coulombic collisions instead of a Foklcer-Planck 
treatment. 
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Fig. 17. Comparison of the effect of using different Coulombic slowing 
treatments . 
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F. Comparison of Methods 

Although both the analytic results of Eq. CIO) and the Monte Carlo 

simulations of Eq. (13) have been presented and common trends observed, 

they have not been compared with each other. In Fig. 18, a series of 

graphs at various temperatures and a constant neutron flux of 2 x 10^^ 

2 

neutrons/ (cm -sec) illustrate the differences between the analytic (dashed 
line) and Monte Carlo (A*s) results. The corresponding tail of the 
Maxwellian distribution (solid line) is included in each of the graphs 
of Fig. 18 in order to locate the intersection of the two distributions; 
i.e., the range of validity of the slowing down distribution. 

The most significant difference in the two sets of results is the 
increasing gap between them with increasing temperature. This dis- 
crepancy can be explained by first observing that if both solutions were 
extrapolated to higher energies, eventually the analytic result would 
intersect the Monte Carlo result and finally lie below it, as is easily 
seen to be the case for 5000“K, From this, it can be concluded that at 
the origin or highest energy for which the nascent source- exists (~2.1 keV), 
the analytic solution will lie closer to but always below the Monte Carlo 
solution, as the temperature is increased. At lower temperatures, the 
importance of the inelastic collisions as an energy loss mechanism in- 
creases, thereby, rendering the assumption of continuous slowing (and the 
analytic treatment) to be invalid at low temperatures. Although the 
inaccuracy may seem insignificant at the origin of the calculation 
(•v2,l keV) , the propagation of the error amplifies the inaccuracy in a 
peculiar manner as described in the ensuing pages. 
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Fig. 18. A comparison of Monte Carlo results (A's) and analytic solution 
(dashed line) for various temperatures at a neutron flux of 
2 X 10^2 neutrons/ (cm^-sec) , Maxwellian distribution is solid 
line. A vertically displaced analytic solution is also presented 
(— ). 
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By examining Eq. (4), the discrepancy in the distribution fvinction 
at the highest energy can be traced to the energy loss rate. Neglecting 
recombination, the equation can be rewritten at the point of origin of 
the calculation subject to the boundary condition 


CE-H-dLE.) = G 

as. 


C39j 


const oLn-b = ScEnAe = -fcE:) ^1 

At 


(40) 


Eq. (40) demonstrates that an exaggeration of the energy loss rate 
will result in the underestimation of the distribution function as is ob- 
served to be the case with the analytic result in comparison to the Monte 
Carlo result. 


The reason that the analytic method overestimates the energy loss 
rate lies in the continuous slowing approximation as demonstrated by the 
following examples. 

Consider an energy cell of width AB located at the point of origin 
E + dE of the calculation, containing 100 particles distributed evenly in 
energy, i.e., 10 particles are in the subinterval labeled A, etc. Csee 
Fig. 19). Furthermore, assume that an arbitrary fraction of the total 
particles collide per At, e.g., 1/10. If the average energy loss per col- 
lision is less than AE, e.g., AE/2, then the number of particles which 
leaves the cell due to collisions according to both treatments is: 
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where represents the. one particle of the ten particles in sufainterval 

A that collides, losing sufficient energy to escape the interval AE within 
the time At. Hence, both treatments yield the same result if the energy 
transfer is indeed infinitesimal. However, if the average energy lost 
per collision is larger than AE, e.g., 2AE, then 
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and 
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The results of Eqs. (45) and (46) confirm that the analytic result errone- 
ously overestimates the diffusion from the original cell by spreading out 
the energy loss per collision over all' of the electrons within the cekl, 
enabling more of them to leave. 

At inteimiediate energies, both the number diffusing into an energy 
interval dE about E and those diffusing from the interval will be erroneously 
calculated in the analytic method. Because the cross sections are rela- 
tively constant and the energy loss per collision small compared to the 
electron energy, the errors will cancel, yielding an approximately correct 
slope for the distribution function at intermediate energies. This is bom 
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out by the result at 5000“K and a renormalized result fdot-dash line) 
at 8000“K in Fig. 18. 

At lower energies, i.e., near the point o£ intersection with the 
Maxwellian, the assumption of continuous slowing down completely breaks 
down. The inelastic cross-sections vary rapidly and the energy lost 
per excitation collision becomes a sizable fraction of the electron energy. 
Thus, the shapes of the Monte Carlo curves differ considerably from the 
analytic results in the region at, and below, the intersection with the 
Maxwellian. Fortunately, the worse departure occurs below the inter- 
section where the calculation is no longer interesting. 
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CHAPTER IV 
CONCLUSIONS 

A. Review 

Preliminary to the’ calculation of the distribution function in a 
uranium plasma is the assemblage of the excitation and ionization cross 
sections for neutral and singly ionized uranium. As no previous calcu- 
lations nor measurements exist for them, the hydrogenic model of 
Gryzinski is applied to uranium utilizing the atomic state data of 

ri71 

Parks. •' Furthennore, the nascent electron source, heretofore an 
undetermined quantity, is modeled upon a semi-empirical formulation of 
the fission-fragment thermalization process and the aforementioned cross- 
section set. Then, these calculations in conjunction with the determination 
of the species' densities via the Saha equations serve as the basis for 
the distribution function calculation. 

The distribution is decomposed into two parts: a Maxwellian (valid 

at low energies) and a high-energy tail. The calculation of the tail is 
performed via two distinct methods. The first method was based upon the 
assumption of continuous slowing down and yielded an analytic solution 
from which trends could easily be predicted. The second method involved 
the Monte Carlo simulation of a governing equation in which the assumption 
of continuous slowing down had been relaxed for inelastic collisions. The 
second method affords a check upon the first method while serving as a 
powerful tool for performing detailed calculations. The disparity in the 
two sets of results is tracable directly to the applicability of the con- 
tinuous slowing-down assumption to the inelastic collisions. 
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Incorporated into the Monte Carlo simulation are several noteworthy- 
techniques developed or adapted specifically for the present case. Fore- 
most of. these are the adaptation of the rejection technique to increase 
the sampling efficiency and the development of the convergence propagation 
tefchnique and the scheme of superimposing continuous and discrete slowing. 
Especially important is the latter technique as it represents the first 
time that Coulombic collisions are considered in irradiated plasmas which 
are being examined for their excitation capabilities. 

In an attempt to ascertain the effect of collective interactions 
upon the Coulombic energy loss rate, two theories were employed. The 
first, the Fokker-Planck theory, is a zero order treatment in the plasma 
parameter g. The other, the unified theory, is a first order theory. On 
the basis of these two theories, it was concluded that the collective 
interactions were adequately incorporated into the calculation. 

From the results presented here, the distribution function is con- 
cluded to be non-Maxwellian above 15 eV. Parametric studies reveal the 
amplitude of the high-energy tail to be linearly proportional to the 
neutron flux level and inversely proportional to the temperature. The 
degree of deviation of the high-energy tail from a Maxwellian, can be 
gauged by -the following example: for the plasma conditions of 8000“K and 

2 X 10^^ neutrons/ (cm^- sec), the calculated distribution induces 6 x 10^^ - 
more excitation events/ (cm^-sec) than a Maxwellian distribution. 

B, Accuracy of Results 

From Section F of Chapter III. it can be concluded that the Monte 
Carlo solution gives a much more accurate account of the collisional 
processes than does the analytic solution, hence, a more realistic dis- 
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tribution function. The Monte Carlo result can be considered statisticallx 
uncertain to within fifteen percent based upon the results of Table 3. 
In^iroved statistics can be affected primarily by increasing the number of 
simulation particles. However, two other factors have a greater influence 
upon the ultimate accuracy of the results; namely, the errors inherent in 
the cross-section set and the Coulombic collision treatment. 

The other researchers have successfully applied 'the hydrogenic model 
of Gryzinski to calculations of cross sections (e.g., Lo^^^^ employed, 
the' model for Helium) . However, comparisons of the calculated cross section 
to cesium, which is similar in electronic structure, reveal that the calcu- 
lated cross-section set may be somewhat low. The uncertainty involved is 
estimated to be a factor of 2. 

The uncertainty associated with the Coulombic energy loss rate is 
difficult to predict. The error estimating scheme of the BBGKY hierarchy 
is not applicable as the plasma parameter g is not small (i.e., g £ 5.2). 
However, the ratio of the unified slowing expression -Cexact to order g) 
to the Fokker-Planck expression (accurate to order 1) is found to be only 
of order 3/2 in spite of the two treatments* diversity. The Fokker-Planck 
expression depicts a test particle's interactions with the background, 
within the annulus defined by b < r < X.., where b = and X„ are the close 
impact parameter and the Debye length, respectively; while the unified 
theory depicts those interactions within the annulus O' < r < ». Because 
each treatment is so different, yet their results are in good agreement, 
it seems reasonable to conclude that the actual energy loss rate is close 
to that predicted by these two theories. If additional correlation (or 
collective) effects enter, the energy loss rate would be even larger than 
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that predicted by the unified theory. A liberal estimate of the energy 
loss rate would be a factor of 5 larger than the Fokker-Planck result. 

The uncertainty associated with the inelastic cross-sections can 
increase the distribution by a factor up to two particularly at high T. 
Similarly, the Coulombic energy loss rate could possibly introduce a fac- 
tor of two in the opposite direction also preferentially at High T. Then, 
the effects tend to cancel, and the final result can be more accurate 
than either of the components of the calculation, i.e., the distribution 
functions reported here Pigs. 14 and 15) are uncertain to within a 

factor estimated to be considerably less than two. 

C. Implications to Uranium Plasma Program 

The results presented in Section E of Chapter II and III clearly 
demonstrate the distribution function to be non-Maxwell ian. The importance 
of this is partially lost since the bulk of the excitation out of the 
ground state is done by electrons below the intersection in energy of the 

tail and the Maxwellian. For. 8000°K, a Maxwellian distribution will cause 

24 . 3 I 

approximately 1.4 x 10 excitations/ (cm -sec). Above 22 eV, however, the 

13 3 

Maxwellian will cause 3.2 x 10 excitations/ (cm -sec) while the high-energy 

14 13 3 

tail will cause 6.6 x 10 and 3.9 x 10 excitations/(cm'^-sec) for neutron 
fluxes of 2 X 10^^ and 2 x 10^"^ neutrons/ (cra^ -sec), respectively. These 
excitation rates may not produce inversions of the excited state densities 
in a uranium plasma. However, if the plasma were seeded with a species with 
a high threshold energy for the first excited state, e.g., helium, then 
ideal conditions exist for predominately exciting the helium with the non- 
Maxwellian tail. Then, clearly, the significance of the high-energy tail 
is determined by the type of interaction under consideration. 
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D. Future Work 

1. Further Applications and Development of Code 

The presence dependence of the distribution should be investigated. 
Crucial to this investigation is the size of the perturbation to the Saha 
densities caused by the fission-fragments. The changes in the densities 
of the various plasma species will determine if the distribution can be 
extrapolated from the results presented here. If not, the calculations 
must be repeated with the appropriate densities (and corresponding tem- 
peratures) . 

' An additional effect that should be incorporated is spatial diffusion. 
This effect is neglected in this work as the projected size (1 meter) of the 
device containing the uranium plasma makes the plasma an infinite one. 

The electrons in the high-energy tail will- preferentially leak out from 
the system, thereby decreasing the perturbation of the high-energy tail on 
the Maxwellian. 

One must not neglect the need for implementing improved inelastic 
cross sections as they become available. Also, the Coulombic inter- 
actions are in need of a model which can more precisely describe the 
interactions under the unique plasma parameters incurred. By far, these 
two aspects are most in need of improvement. 

2. Analysis of Other Plasmas 

From the conclusions reached in section B, the code needs to be- 
expanded in order to accommodate the presence of seed species as well as 
buffer gases. These additional species introduce important factors in 
the calculation of a distribution function in a working uranium plasma 
and the initially planned experiments. Their presence will dilute the 
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•electron source and in the case of the seed gas decrease the stopping 
power of the plasma. It is difficult to judge a priori if these opposing 
trends will cancel each other. Furthermore, the introduction of molecules, 
e.g., UFg into the plasma will complicate the calculation of energy loss 
rates -with the introduction of vibrational and rotational- excitation 
processes. The ease with which these excitation processes occur will 
increase the energy loss rate, decreasing the high-energy tail. A more 
detailed analysis is possible through the adaptation- of existing treat- 
ments of electron slowing via molecular excitation collision. 

The code might also be applied to other plasmas with a distributed 
source of nascent electrons. The plasmas will have to be restricted to 
those which satisfy the assumption that the high-energy tail is only a 
perturbation to a Maxwellian distribution. This excludes plasmas with 
electric fields wherein the bulk of the plasma is describable by a 
Druyvestyen distribution. 
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APPENDIX A 
ENERGY LOSS RATES 

1. Fokker-Planck 

For the majority of the calculations reported in Chapters II and III, 
a traditional Fokker-Planck expression for the Coulorabic energy loss rate 
was employed. The limited applicability of the Fokker-Planck expression 
due to the small number of electrons within a Debye sphere was pointed 
out previously in Chapter I. However, the expression was used in spite of 
the expansion factor g not being negligible compared to one, following 
the lead of others in applying the expression under similar circumstances, 
e.g., in MHD calculations of conductivity and in modeling afterglows . 
The precise equation utilized is; 




(47) 

where n^, m^, q^, and T^ denote the background species s's density, mass, 
charge, and temperature; ra, q, and v denote the test particle's mass. 


charge, and velocity; and the Debye length is The function F 

D s 

appearing in Eq. C47) is defined as: 


Calculations of the energy loss rate as prescribed by Eq. (47) appear in 
Figs. 20 and 21. A discussion of the results appears in the next section. 
2 , Unified Theory 

In Chapter I, the speculation that -collective interactions play a 
dominant role in the slowing of energetic electrons leads to the need to 
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consider a higher order kinetic equation, i.e., the Lenard-Balescu^^^^ 
equation. The traditional treatments of charged particle slowing in the 
Lenard- Bales cu ox wave theory, e.g., Sigmar and Joyce, introduce a 
cutoff to eliminate a singularity in the analysis. The uncertainty as- 
sociated with the introduction of a cutoff, especially under the plasma 
conditions present in the uranium plasma, render these treatments inap- 
propriate. A theory which is independent of a cutoff was developed by 
Kihara and Aono.^^^^ Since this theory depicts the full spectrum of 
interactions, i.e., from collective interactions to close impact col- 
lisions, it is exact to order g. Then, through its implementation, the 
effect of neglecting collective interactions can be gauged. 

The theory of Kiharo and Aono or the unified theory is based upon 
the observations of Hubbard that the divergences appearing in traditional 
wave and impact theories could be made to cancel each other. Symbolically, 
this can be written as 


Where X denotes a reaction rate, e.g., the n^^ moment , in either 

the close impact or wave theory. The subtrahend 2ui expression which 

neglects the effect of collective interactions and the effect of orbital 


curvatures, is responsible for the cancellation of the divergence in both 
these theories. Although Hubbard's formulation also includes the full 
spectrum of interactions, the resulting energy loss rate is not indepen- 
dent of cutoffs. 


The prescription for calculating a relaxation rate in the unified 
theory is dependent upon casting a relaxation rate in both the wave theory 
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and impact theory in the following forms 


X = Kck.) 


(50) 




and 


(51) 


3<Ab Bcb) 

b>0 

where k is the wave number in the wave theory and b is the impact parameter 
in the impact theory. The theory of unification then states that the non- 
diverging relaxation rate is given by 

O o 

where b^ is an intermediate length, less than the Debye length and greater 
than the close impact radius. A relaxation rate calculated according to 
Eq. C52) is independent of the intermediate length b^. 

As an exercise in the unified theory, Itikawa and Aono^^^^ calculated 
the relaxation of a test particle of arbitrary velocity in a plasma. Their 
result is of the foonn 


dt 






(53) 


where Any = 0.5772, v^ = 



s ”^ 5 "^ 

> ““ - -STTST 


The term G contains 
s 


s '•■■'s 

the Coulomb logarithm dependence upon the velocity of the test particle. 

The exact form of G is^^^^ 

s 

G,;(u)= Q(U) - Cjj^cu) C54) 
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where 


u = v/v and 
s s 


>0^ 




QCu) = a Qcu) +- 


4>cu)-- u 


z,u 


C55D 


Gcu) = 


4>(u)~ gcjjCu) 


2- u' 


C56) 


Cu) — ^ 6 




(57) 
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(58) 


=vt j e C,,m 


(59) 


-U‘ 
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(60) 


C ,(.u) = 

3S^ 


^ Jn.( /)„.(u) + B,^,Cu)) 


+ 


2| B,,(u)l 


(f- 


I B„,(u)l , 


(61) 


A..<-r) = + As- ^ 


(62) 
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c< r= Vs 

/T.m, 
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Due to the complexity of Eq. (54) in its entirety, a more tractable 

f391 

expression for G^, as developed by Perkins, ^ was substituted in place 
of the exact expression. Perkins* expression [Eq. (68)] is based upon the 
asymptotic behavior of Eq. (54) i.e.. 




where 


= 3i^ — '/z. 

= JnZ - c 


= 2.in-^ v>v, 

= X>v>v, (70) 


and the subscript 1 denotes the electrons and the 2 denotes the ions. 
Some error is introduced into the energy loss rate when Eq. (68) is em* 
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ployed as v approaches V 2 . However, this error can be gauged to be 
negligible for the electron energies considered through an analogous 
•con^iarison of a Butler and Buckingham type expansion to the Fokker- 
Planck expression. 

The Coulombic energy loss rates predicted by both the Fokker-Planck 
theory and the unified theory appear in Figs. (20) ^d (21) for plasma 
temperatures of 5000“K and 8000'’K respectively. The disparity between 
the two theories grows as the electron energy decreases until a maximum 
is reached between 2 and 3 eV. This effect can be attributed largely' 
to the better coupling of the electrons with plasma waves as the electrons 
approach the wave velocity. By increasing the temperature, the disparity 
is observed to increase only slightly. This is the result of the increased 
presence of collective interactions due to a further penetration into the 
classical collective plasma region of Fig. 2. 

The history of a test electron as it slows down in each of the slow- 
ing, theories is displayed in Figs. 22 and 23. The disparity in the slow- 
ing profiles is attributed to the increased stopping power in the unified 
theory due to the presence of collective interactions. A measure of this 
disparity is the thermalization time, defined here to be the time to slow 
from 1 keV .to 1 eV. From Figs. 22 and 23, the Fokker-Planck theory is 
observed to yield a thermalization time approximately 1.6 times larger 
than that predicted by thetmified theory. The factor 1.6 represents 
a first approximation to the inaccuracy introduced by the small number' of 
particles within a Debye sphere. Although this degree of inaccuracy would ' 
seem undesirable, the inaccuracy associated with the inelastic collisions, 
which are equally important in determining the total energy loss rate, is 
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ENERGY, eV 

Fig. 20. The Coulombic energy loss rates for an electron slowing off of 
the various plasma species and the total loss rate versus energy 
for both Fokker-Planck (FP) and Unified Theories (UT) in a 5000“K 
plasma. The total and electron energy loss rates are denoted by 

C )FP and (-) UT, The singly ionized uranium energy loss rate 

is denoted by ( )FP and ( )UT, while the doubly ionized 

uranium energy loss rate is given by ( )FP and C-**-!)UT. 
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ENERGY, eV 

Fig, 21. The Coulombic energy loss rates for an electron slowing off of 
the various plasma species- and the total loss rate versus , 
energy for both Fokker-Planck (FP) and Unified Theories (UT) 
in a 8000*K plasma. The total 'and electron energy loss rates 

are denoted by ( -)FP and C-)UT* The singly ionized uranium 

energy loss rate is denoted by _C )FP and C )UT, while the 

doubly ionized- uranium energy loss rate is given by (••••)FP and 
C 3UT. 
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TIME, sec 

Fig. 22. The energy of a test electron is plotted versus' time as it slows 
according to FP and UT theories in a 5000® plasma -denoted by 
( — ) and C-i*.-) respectively. The energy gained by each of 
the plasma species i‘s also plotted versus time with the same 
delineation of species and theories as in Fig. (20) , 
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TIME, sec 

Fig. 23. The energy of a test electron is plotted’ versus time as it slows 
according to HP and UT theories in a 8000*K plasma denoted by 
C — ) and respectively. The energy gained by each of 

the plasma species is also plotted versus time with the same 
delineation of species and theories as in Fig. C21). 
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even larger. Then, for the present calculation, the agreement between 
the unified theory and the Fokker-Planck theory is satisfactory 
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APPENDIX B 
CROSS SECTIONS 

Appearing in Fig. 24, are the ionization and excitation cross, sections 

i 

for electron’ bombardment of neutral, singly, and doubly ionized uranitim. 

The cross sections have been calculated from formulae based upon a sym- 
metrized version of the Gryzinski model as reported by Burgess and 
Percival, implementing the ionization and excitation data of ParKs, 
et>al.^ Because of the lack of a well-defined state, corresponding to 
the ionization potential of the various uranium species, it has been 
assumed that a limited number of the outermost electrons-, usually 8 or 
less, participate in both ionization and; excitation processes. The ex^ 
citation cross section represents a sum of cross sections for transitipns 
from the ill-defined ground state, consisting of any of the outermost 
electrons employed in the representation of a state corresponding to the 
ionization potential, to a multitude of excited states, 27 total excited 
states for each species (see Table 4) . The transitions are governed by 
the selection rule |A2.[ = 1. The multiplicity of the allowable transi- 
tions eliminates the resonance behavior typically exhibited in excit.ation 
cross sections, e.g., cesium. 

As can be seen from Fig. 24, the cross sections exhibit an abrupt 
rise at the threshold energy as is characteristic at the onset of a 
quantum mechanical process. For energies below the threshold, the cross 
sections approach zero, -but because it is not of immediate interest, this 
threshold region has not been explored in greater depth. The appearance 
of discontinuities in the slopes of the cross-sections are indicative of 
one electronic state's participation in an event being overshadowed by 
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Table 4. The uraiiium states and their corresponding quantum numbers 

employed in the cross-section calculation (after Parks, et al 
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Fig. 24. The ionization and excitation cross sections for U, U' 
U'*"'' versus energy. 
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another state; e,g,, the discontinuity in the slope of the ionization 

cross section at 18 eV represents the' appearance of an inner-electron 

ionization process which, at high energies, overshadows the ionization of 

the outer-most electron. At very large energies, the fine structure of 

the atom gives way to an 1/E energy dependence. 

The excitation cross sections are found to* be larger in magnitude 

than the corresponding ionization cross-sections. A. similar trend is ob- 
041") 

served for .cesium^ ^ which is quite similar in electronic structure. 

These observations can be attributed to a cross section, in general, being 
inversely proportional to the energy transferred. Then, the smaller- energy 
transfer afforded by excitation collisions result in the increased probabil- 
ity of their occurance over ionization collisions. 

A comparison of the magnitudes of the uranium cross sections presented 
here to cesivmi ionization and excitation cross sections reveals that 
the uranium cross sections are smaller by an order of magnitude than the 
cesium cross sections which are measured accurate to within a factor of 2. 
This trend can.be explained by examining their ionization potentials; 
namely, 3.89 eV as opposed to 6.22 eV for neutral uranium. The 

similar electronic structure of the two elements as well as the relatively 
close proximity of their ionization potentials suggest that the uranium 
cross sections should be somewhat smaller due to uranium’s larger ioniza- 
tion potential (implying a larger energy transfer) . This is in qualitative 

agreement with the present results. However, a comparison with helium 
(42) 

cross sections reveals the uranium cross sections to differ from the 
cesium cross sections by a wider margin than anticipated. The uranium 
cross sections are an order of magnitude lower than cross sections of 
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cesium and other elements with similar ionization potentials. .Yet, the 

Gryzinski model was found to be inaccurate fay at most a factor of two for ‘ 
f22'l 

helium. ^ Then, without further ejqserimental data for guidance, a crude 
estimate. of the error introduced by applying a hydrogenic model for cross, 
sections to uranium can be ascribed to be a factor of two too low. 

The ability of an ionized particle to focus approaching electrons, 
thereby, enhancing the cross section is also taken into .account in these 
calculations in the manner prescribed by Burgess and Percival. This effect 
can best be seen by comparing the excitation cross section of U and U^. 

One would expect the tighter bound electrons of to be harder to ex- 

cite. However, the charge on U not only draws in electrons that would 
normally pass on by, but it also gives them additional energy as they are 
accelerated in the potential field. This more than compensates for the 
tighter shell structure of U^. 
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APPENDIX C 
COMPUTER CODE 

A digital computer code was developed to solve Eq. (26) employing 
the simulation technique described in Section C of Chapter III. The 
simulation technique was augmented by various variance reduction tech- 
niques to enhance the rate of convergence of the solution. In order to 
further reduce conqiutational time, several computational "tricks'' were 
employed. All of which are to be described in the ensuing pages. 

The computer code can best be described through references to the 
flow chart in Fig. 25. Although the flow chart provides an over-simpli- 
fied view of the program, the spirit of the calculation is preserved by 
it. Nvimerals have been placed adjacent to the flow diagram to aid in 
the identification, of various sections of the code. Let us begin to 
follow the flow of the logic with Section 
Input 

The first section contains the input parameters. One such parameter 
is the random number starter. By altering this niimber, a different ran- 
dom number string is used as the basis of the Monte Carlo simulation. 
-Such freedom is essential to statistical testing- of the results. 

The plasma properties such as density, ten^jerature, neutron flux, 
and identification of background species through mass and charge con- 
stitute a second set of input parameters. These parameters permit para- 
metric studies of the effect of the plasma environment upon the relaxa- 
tion of nascent electrons into a Maxwellian distribution. Some of the 
parameters are correlated, such as temperature and density. These para-- 
meters must be self-consistent upon input as they detemnine into which 
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Fig. 25. Flow Chart of Monte Carlo Code 
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specific target Maxwellian the nascent electrons will relax. Although 
the plasjna species can be changed through the input data, the change is 
not complete^ without providing a new subroutine containing the appropriate 
cross --section set. 

Another input parameter consists of an initial guess for the dis- 
tribution function. Previous attempts at Monte Carlo simulations of dis- 
tribution functions conclude this to be a key element in reducing compu- 
(231 

tational time. The present program was not very sensitive to this 

feature. This can be attributed^ to the propagation of convergence tech- 
niquC' whose application is made possible by the near linearity of the 
problem. 

Also, the nascent electron distribution function is required as an 
input parameter. The nascent electrons constitute only a portion of the 
total electron source. The remaining source, namely the secondary elec- 
trons, is calculated within the program and is consistent with the col- 
lision rate. As before, the nascent source must be self-consistent with 
the other input parameters. 

Finally, the convergence criterion and the number of desired itera- 
tions of the converged solution must be supplied as input parameters. 

The convergence criterion in conjunction with the number of simulation 
particles determine the accuracy of the solution and will be explained in 
more detail later. The parameter for the number of iterations provides 
a means of data smoothing. It essentially determines the number of 
snapshots of the distribution in time which are to be used in formulating 
the basis of an average. An averaging process is necessary to rid the re- 
sults of fluctuations which are characteristic of all Monte Carlo Simula- 
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tions (see Figs. 11 and 13). 

, 2. Initialization 

a. Fitting Functions of Frequently, Calculated Expressions 
The second section of the program contains those operations termed 

initialization. These operations can be divided into two types, time- 
saving and preliminary calculations. The distinction between them is 
that the preliminary calculations must precede the remaining sections 
of the code, whereas, the time-saving calculations are more conveniently 
calculated at the earlier stages of the program so that duplicate calcu- 
lations may be avoided. The tabulation of integrals and the cubic-spline 
fitting of frequently needed, complex functions (i.e., the expression for 
E(t) plotted in Fig. (10)) are examples of time-saving calculations, 

b. Allocation of Machine Particles 

The preliminary calculations entail the distribution of the machine 
-particles amongst the energy regions for which values of the distribution 
function will be calculated. The choices of average machine particle 
density and density distribution affect the precision of the calculated 
electron distribution function both globally and locally. The ideal man- 
ner in which to distribute the machine particles would be to mimic the 
expected electron distribution. However, the range of the variation in 
electron density and the fact that calculation of the extreme lower end 
of the distribution is unnecessary make such an approach impractical. 
Since the motivation of this work is to provide a basis for calculating 
excited state densities, the emphasis should be placed upon excitation 
rates, which suggests, the distribution should mimic the total macroscopic 
excitation cross section, i.e., the range of energies where the highest 
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degree of accuracy is obtained contains the most particles which con- 
tribute to the excitation rate. This insures a high degree of accuracy 
in any calculation of excitation rates based upon the results for the 
distribution function obtained from this code. Due to the wide range of 
values of the cross section over the energy range of the calculation, 
these variations needed to be toned down. The distribution finally ar- 
rived at is 

M(0=c+i^XlCE) ™ 

where c is a constant determined from the average error permitted and 
2(E), the total macroscopic cross section at energy E, approximates the 
excitation macroscopic cross section. For 1700 machine particles, Eq. (71) 
yields a minimum of 50 particles per group and a maximum of 135 particles 
per group. 

3. Convergence Propagation 

The third section of the code denotes the implementation of the 
convergence propagation technique. The method is to take advantage of 
the dependence of f(E*) upon only f(E>E'), since the interaction of 
f(E*) with f(E<E') is negligible compared to its interaction with the 
Maxwellian part of the distribution. Hence, it is inefficient to simu- 
late the lower end of the tail of th? distribution wh^e simultaneously 
simulating the upper end of the tail, if the results for the upper end 
have not yet converged- to the final solution. Thus, the calculation be- 
gins at the upper energy region until convergence is obtainedj and then 
the simulation is expanded lower in energy, one region at a time. Such 
an expansion in energy is analogous to the propagation of a wave, from 
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one group to another, whence the name. 

Convergence of a group is established by matching the flow of par- 
ticles into a group with the outward flow from the group within a specified 
tolerance level. Care must be taken not to set the tolerance level be- 
low the noise level of the Monte Carlo simulation. The noise level is 
a manifestation of the fluctuations characteristic of all Monte Carlo 
simulations and can be approximated by the square root of the machine 
particle density within an energy group. If precautions are not taken 
and the tolerance level is set below the noise level, a superficial con- 
vergence is obtained through the compounding of random fluctuations in- 
particle flow, producing an vtnphysical result. 

4. Secondary Initialization 
a. Initial Distribution 

Appearing in the fourth section of the program is a belated 
initialization phase. Here, the cujnrent lowest energy group is initial- 
ized to a guess distribution. In so doing, the computational time re- 
quired to obtain convergence is minimized to the degree to which the guess 
approximates the solution. Computational time is further reduced by in- 
troducing into an energy group k • m particles at m discrete energies and 
randomly staggering their associated times. This permits the repetitive 
use (k times) of the various probabilities CHqs. (35) and (37)) necessary 
for the determination of the velocity of the next collision and the par-- 
tide weights. The assignment of particle weights is straightforward, i.e., 
the number of electrons introduced into an energy interval, determined by 
the guess distribution, are evenly distributed over the machine particles 
assigned to the interval. 
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b. Varying Time Step At 

The time period for which the particles are permitted to evolve is 
also determined in this section, permitting the period to be varied as 
dictated by convergence and efficiency requirements. Since the time 
width of an energy group (or the time for an electron to traverse an 
energy group) is progressively smaller for decreasing energies, the time 
period can be determined solely by the lowest energy group involved at 
any moment during the calculation. The upperbound on the time period must 
be less than the time width of this group in order to fully utilize the 
initial guess. A lowerbound is established by the graininess of the 
guess (or the number of discrete energies m at which the guess particles 
are stacked). For intermediate values, the time period is determined by 
requiring the number of source particles introduced into the group to be 
a fraction of the total number of particles within the group. Hence, we 

have . ^ see.-) ^ 

' -Pee) 

-fcE) P2) 

■5'SCE-)'f 

where T is the time width of the current lowest energy group, f(E) is 
the initial guess for this group, S(E) is the nascent source rate, and 
At is the calculated time period (the same At as in Fig. 10). The 
nascent source rate is used to approximate the total source rate in 
Eq. (72) to speed up the calculation. This is possible since only At, 
not the calculation itself, is affected. 
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5. Processing 

a. Particle Update 

In the fifth section of the program, the actual processing of par- 
ticles is performed. Initially, a computational clock is incremented 
by the time step At which is calculated in the previous section of the 
code. Symbolically, this is representative of the next step, namely the 
updating process wherein all particles will be permitted to evolve for the 
period At. The evolution of the individual particles is affected as pfe,- 
viously outlined in Section C of Chapter III, where the particle Simula- ■ 
tion is described in detail. The simulation algorithm is programmed as 
outlined in the flow chart in Fig. 26. 

The flow chart indicates that the coding requirements to handle each 
of the three collisional types and hence, the costs of computation, gets 
progressively larger when proceeding from Goulombic interactions to ex- 
citation collisions aiid even larger in going to ionization events. An 
explanation for this trend follows. The inelastic collisions require ad- 
ditional coding to arrive at the electronic state in the target that will 
participate in the collision. The energy lost due ,to.^ excitation col- 
lision. can be estimated to be the energy of excitation, thereby minimizing 
costs. For an ionization collision, the energy lost must be calculated by 
a rejection technique (subroutine EP.S) . Furthermore, supplemental coding • 
is necessary in the advent of the birth of a secondary electron energetic 
enough to influence the calculation. This portion of the code is dominated 
by the time spent in calculating the energy where the secondary electron 
first. collides. A similar calculation of the next collision energy must 
also be performed after the occurrence of either type of inelastic col- 
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Fig. 26, Enlargement of the section of coding which updates particles, 
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lision. This calculation is one of 'the most frequently and most costly 
performed operations in the Monte Carlo program. The algorithm employed 
is a form of the rejection technique. The logic of the technique was 
described in Section C of Chapter III. Due to the possibility that a 
candidate for the collision energy may be rejected, the algorithm is a 
potential infinite loop (see Fig. 27) . The number of circulations through 
the loop can be minimized by retaining one of the random number pair con- 
nected with the rejected candidate. Specifically, the random number 
is compared in the rejection decision process to the probability of a 
collision at the candidate energy determined by r^. Then T 2 is subtracted 
from one and used as r^^ in the next loop- The process can be viewed in 
the dart analogy (see Section C of Chapter III) as reflecting the point 
of impact of a rejected dart from the upper right hand comer to the low- 
er left hand comer where the dart ivill more likely fall below the prob- 
ability distribution, resulting in the acceptance of the candidate for 
collision energy. The mapping of T 2 onto r^^ does not generate a true 
reflection. However, the collision frequencies were not altered by the 
above scheme. A similar scheme was also employed in the rejection tech-, 
nique used in calculating the energy lost due to ionization collisions 
with the same degree of success. Furthermore, the scheme proved to be 
more efficient than that used by Carter, et al . 

During, and after the particle has been updated, several observa- 
tions are recorded pertaining to its activities. These include the 
initial and final location in energy space, the energy lost as well as 
the type of collision responsible, and the production (if any) of secon- 
dary electrons. These observations serve as the raw data for the conver- 
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Fig. 27. Flow chart of the rejection technique employed in calculating 
the energy of the next collision. 
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gence test and the final results. 

b. Secondaries 

The registering of the production of a secondary electron entails 
the creation of a new particle with the sa^e weight ‘as its parent. Like 
its parent, it must' complete its evolutionary period At from the time of 
the collision, i.e., the time At* in Fig. 10. By offsetting the time of 
. each secondary at birth by the quantity At '-At, the entire batch of 
secondaries can be updated as a group for time At after the original 

are processed 'with no distinction made in the updating process 
as to the particles origin. Similarly, succeeding generations of secon- 
daries can be generated and processed, culminating in the cessation of 
the avalanche of secondaries. 

c. Source 

i. Vacancies 

After secondaries are introduced into the electron population, the 
nascent electrons are generated as source particles. Since the machine 
particle population in an energy group is fixed, the injection of source 
particles into the particle population requires two algorithms to accom- 
modate both vacancies and an excess of electrons in each group. Vacancies 
are filled with the nascent electrons in the same manner that a vacant 
group is initially filled with a guess distribution. However, the par- 
tide stacking procedure is complicated by the fact that the ntimber of 
vacancies is not always factorable into the product of two integers k and 
m, where k is the number of particles to be stacked at the ra discrete 
energies within the energy group, as was previously done. Nevertheless, 
the number of vacancies is resolvable as k • m + A, where A is also an in- 
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teger, thereby permitting a partial realisation of the advantages of 
stacking particles. 

ii. Excess 

Should there be an excess of electrons in a group, a source particle 
must displace the excess particles in order to be accommodated into the 
group Csee Fig. 28). Immediately, the question arises as to which par- 
ticles in the group constitute the excess destined for extinction. This 
can be resolved by seeking those particles whose removal will yield the 
least perturbation on the calculation, i.e., those with the least weight. 
Once assembled, they are eliminated while retaining knowledge of the 
weight of the ensemble and the fraction .which were to have collided. 

This information is combined with the weight of the source particle to 
be introduced and the corresponding probability of incurring a collision 
to generate a hybrid source particle to be added to the energy group. In 
this manner, none of the information retained by the "killed" particles 
is lost. 

6. Output 

In the sixth and final section of the program, the individual 
particle observations previously defined are accumulated for the calcula- 
tion of the following quantities: flux into and out of the currently 

lowest energy group, secondary electron production rates, energy loss 
rates, the W-value (defined to be the energy lost by electrons per ion 
pair formed), and, of course, the distribution function-corresponding ‘ 
to the time displayed upon the internal clock for the energy range for 
which convergence has been established. The particle fluxes and secondary 
production rates are combined with the nascent source rate to yield the 
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Fig, 28, Enlargement of the section of coding which kills off the excess 
particles in an energy group. ■ 
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particle balance in energy space which serves as the convergence criterion. 
The particle balance, energy loss rate,, and' W-value .are secondary output, 
displayed as a check upon the validity of the results and for future use.' 
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